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I. NTRODUCTICH 


The investigation of wave rotors! and their possible 
application in gas turbine engines at the Naval Postgraduate 
School's (NPS) Turbopropulsion Laboratory prompted the 
development of a one-dimensional unsteady flow computer 
code. The QAZ1D method developed by Verhoff [Ref. 1] was 
chosen over other methods [Ref. 2] for reasons which 
included the following: 


1. The method is based on the use of characteristics. 
Such methods can model wave propagation accurately. 


2. The use of a natural streamline coordinate system 
eases the difficult task of computing with two and 
three dimensional grids. 


3. The equations are written in a form which allows a 
straightforward extension to viscous flows. 


Salacka [Ref. 2] verified the development of the equations 


and implemented a one-dimensional Euler solution of the 


lfhe wave rotor consists of simple tube-like passages 
along and around the periphery of a drum which rotates 
between end walls containing inlet and outlet (partial 
admission) gas ports. Compression and expansion processes 
are arranged to occur in a cyclically periodic unsteady 
process within the rotor such that a high pressure driver 
gas is used to compress a low pressure driven gas. The 
compressed driven gas is led from the outlet port to the 
combustor and brought back into the rotor as the driver. 
The expanded driver gas exits an outlet port. Thus the 
rotor functions as both turbine and compressor with direct 
gas-gas energy exchange through unsteady wave processes. 
The technology of such devices requires the design of 
appropriate cycles using one-dimensional calculations, and 
analysis of the multi-dimensional unsteady flows such as 
occur during port opening and closing. 
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shock tube problem in the EULER1 code. The EULER1 code was 
limited to high pressure set on the left, with flow and 
shock velocities to the right. The shock wave was tracked 
and corrected for in the flow, but the contact Peconemnni ty 
was not, and the expansion wave was not tracked. The code 
also did not treat end boundary conditions. 

The present work extended the EULER1 code of Salacka to 
become the "Euler 1D Version 2" or E1DV2 code which: 


1) can handle both right and left traveling flow and 
discontinuities 


2) implements tracking and correction of the contact 
discontinuity 


3) implements tracking of the expansion wave 
4) models closed wall and open end boundary conditions 


5) models shock and contact surfaces within a grid 
interval. 


Section II and Appendix A describe the analysis underlying 
the revisions. The Fortran code, E1DV2, is detailed in 
Section III and Appendix C and graphical results of four 
test cases to demonstrate the new features of the code are 
given in Section IV. A discussion of these results and 
limitations of the present code follow in Section V. 
Conclusions and recommendations are made in Section VI. 
Appendix B contains instructions to operate the code on the 
NPS computer and the FORTRAN listing is included as Appendix 


we 
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II. ANALYTICAL AND COMPUTATIONAL APPROACH 


A... QOAZID DESCRIPTION 
The QAZ1D analysis and numerical solution scheme is an 
explicit, non-conservative method that uses extended Riemann 
variables to model the multi-dimensional Euler equations in 
a natural streamline coordinate system. Table I lists a 
summary of the applicable equations which are derived in 
detail in [Ref. 2]. 
1. The Coordinate System 
Figure 2.1 shows the natural streamline coordinate 
system (s,n,m) relative to a fixed rectangular cartesian 
system (X,y,Z). The right-hand orthogonal system moves in 
curvilinear ‘translation along a streamline. The "s! 
direction is measured in the direction of the flow, tangent 
to the streamline. The "n"™ direction is perpendicular wae 
the s direction in a plane perpendicular to the x-z plane. 
The "m" direction is normal to the yvlane containing thea 
and s directions. The angle ¢ 1s the angle between the s-n 
and x-y planes, and the angle 8 is the angle between the 
velocity vector and the x-z plane in the s-n plane. In the 
one-dimensional problem treated here, the s direction is 
such that velocity to the right is positive, and to the left 


1s negative my (Ref. Viz pee) 
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TABLE I 


SUMMARY OF THE EULER EQUATIONS 





Pefinitions 
Speed of sound A =8,/P/o (i221) 
Modified entropy don 
change qdS = - TT (iva } 
so that 
I = $= ee 2Y - in(P/o! )] (ie) 
R, ray —1) 
Extended Riemann 
variables =q+AS (I.4) 
R= q- AS (eS) 
Euler Equations 
3g Be aie a 2 y (ape 2 
+ qss = - Gas - Sp Rg - S41 
ae i * 30 36 Tc 
Y= asis= + cos 9 <S) (1.6) 
4 . ey ae os Bee 
a (q Ay B i A (S vie Eee ee) 

+ ast + cos 0 3 a) 
3S aS (eh 
tos - =? 

2 
come oo - _ A 3 InP (I.9) 
se «Ct «(Os Yq on 

2 3 : 
a6 4,42 - __A 9 InP (I.10) 
ot dS Yq cos6 om 


21 


streamline 





Figure 2.1 The Natural Coordinate System 


2. Extended Riemann Variables 
The "Riemann Variables" in most publications are 


defined as 


re 


and (29 


2 


R> Gat ozo 


where gq is the velocity magnitude, Ais the speed of sound, 
and yis the ratio of the specific heat at constant pressure 


to the specific heat at constant volume for a particularegae 
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fae. Sp ws). Verhoff modified the definition to obtain 


"Extended Riemann Variables," defined as 


tO 
ll 


Sl cp AVS (I.4) 


R= q-AS (i 25) 


where S is the modified entropy given in Table I [Ref. l1:p. 
2a 
3. The Governing Equations 

The Euler equations in Table I were developed in 
detail (Ref. 2:Appendix Aj from the equations’ for 
conservation of mass, momentum, and energy for a control 
volume. The underlying assumptions comprise inviscid flow, 
negligible body forces, no heat transfer, and a perfect gas. 
These assumptions and the definition of modified entropy 
were used in transforming the equations into a form suitable 


for solution using the method of characteristics. 


B. CHARACTERISTIC SOLUTION 

Metcdtetons (1.6), (1.7), (21.8), (I-.9) and (I.10) are a 
system of slightly coupled quasi-one-dimensional partial 
differential equations (PDE) with eigenvalues qt+tA, q, and q- 
A. These equations can be rewritten along the characteris- 
thes in the time-space domain to become = ordinary 
differential equations (ODE) which are solved by quadrature 


and interpolation. The propagation of each characteristic 


ZS 


curve along its particular trajectory is expressed by the 
respective ODE. [Ref. l:p. 1] 


Each in the system of equations is in the form 





OW ge, ees ae (2525 
If w = w(s,t) it is shown by Salacka [Ref. 2:pp. 20-23] that 
_ ds 2s 
A Siar ( ) 
and 
mae (2.4) 
at 


where A describes the characteristic curve along which w 
changes. 

In this section the solution of the 3D Euler equations 
given in Table I is described. The computer code developed 
in the present work, however, was for the 1D case wherein 
equations (1.9) and (1.10) are not required, and only the 
three remaining equations (1.6), (1.7), and (1.8) are 
considered. 

Equations (1.6) to (I.8) for one-dimensional flow may be 


expressed in matrix form by setting 
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w = R [A] = 0 q-A 0 
S 0 0 qd 
and 
- Gis -4p la -yAl, 
Z = AAS - rl, 
0 


Then equation (2.2) becomes, for the equation set, 


Figure 2.2 shows the characteristic curve in the space- 
time domain between two nodes within the computational mesh. 


Equation (2.2) has a solution 


_ 7 ttdt _ 
jw = - dwt f Z dt (2.5) 
ie 
where §w is the change due to time at a fixed location, and 
Aw is the change due to displacement at a fixed location, 


along the characteristic curve. The value of Aw is 


calculated by interpolation through an iterative scheme 


Zz 


t+ dt| H 


aS/ aE 
at 


Interpolated 
poin 


Figure 2.2 Characteristic Curve Within Computational 
Mesh 


using initial guesses for the characteristic slope, %) tnem 


values at time level t. The line integral is approximated 


by 
t+dt B Se ae = 
Est Zz Z ae 
zdt = dS = f{ yds = +(S.-(s,-AS)) 
J re 5 S ,-AS x \ > Jigen 
= 2z2(As/A) 
= 2(6t) (2.6) 
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The spatial derivatives of q and A in Z are estimated from 
values at A and sj. Finally, Ww 1s calculated, and each 
node is updated to the next time level. 

The Courant-Friedrichs-Lewy (CFL) convergence condition 
requires the domain of dependence of the numerical approxi- 
mation to include the domain of dependence of the differen- 


tial equation (Ref. 4:pp. 336-337]. Figure 2.3 shows that 


Figure 2.3 Computational Grid for the QAZ1D Method 


the computational grid for the QAZ1D method satisfies this 
condition by having interpolated initial data points (points 
1, 2, and 3) in between previous solution nodes (points 4, 
5, 6). The qtA characteristic is estimated from values at 
Pomme, 4, gd Characteristic from values at point 5, and q-A 


characteristic from values at point 6. The use of a maximum 


P74 


time step keeps the domain of dependence of the numerical 
scheme just outside that of the actual equations, aiding 
convergence. Unlike finite difference schemes that are 
conservative and require a numerical diffusion or viscosity 
to handle oscillations or smearing, this method is stable 
without numerical smoothing or artificial dissipation [Ref. 
5:pp. 562-583]. However, the equations do not account for 
the discontinuities across shocks and contact surfaces 
properly. Thus after each time step a one-point correction 
technique is used to correct for shocks and contact 


surfaces. 


Ce DISCONTINULTIES 
1. Expansion/Rarefaction Waves 

The head and tail of the expansion wave travel along 
the characteristic as gradient discontinuities. Flow 
velocity, sveed of sound, pressure, and density are continu- 
ous across a gradient discontinuity but the spatial deriva- 
tives are discontinuous. Entropy remains constant through 
the rarefaction wave. The head of an expansion wave moves 
into undisturbed flow with a speed q+A, while the tail of 
the expansion wave has a velocity q-A. No special treat- 
ment of these waves is required; the characteristics method 
accurately tracks and models their influence. Expansion 
waves are generated when two shocks collide, a diaphragm is 


burst, a contact surface and shock intersect, a shock leaves 
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See rrr a oe pe er ee? —reliee id 


an open boundary, or an open boundary has a pressure lower 
than the pressure inside the tube. [Ref. 3:pp. 13-15] 
2. Shocks 

Shocks are created when a diaphragm is burst, 
compression waves generated at an open boundary coalesce 
into a wavefront with sufficient strength, or when a shock 
and contact surface collide. Pressure, velocity, density, 
and entropy are discontinuous across a shock [Ref. 6:pp. 30- 
31]. The equations of motion (1I.6)-(I.10) are not correct 
for calculating changes through shocks. A method to correct 
for a shock is to use the ratio of the jump in the extended 
Riemann variable across the shock to the speed of sound 
downstream of the shock to find the incoming Mach Number 
relative to the shock (w) [Ref. l:p. 3]. Figure 2.4 shows 


the case of a shock headed to the right. The shock velocity 





me leady Ssteqau 


Figure 2.4 Shock Wave with High Pressure on Left 
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(V;,) is imposed on the flow to produce a stationary shock 
condition. Flow variables are then corrected using normal 
shock relations. Figure 2.4 applies to a condition where 
high pressure is on the left. The shock is moving to the 
right at V, into flow with velocity qg. Table II lists the 
equations developed by Salacka [Ref. 2] for relating the 
extended Riemann variable change to Mach Number. Velocities 
are non-dimensionalized by the speed of sound downstream of 
the shock direction, Ap, and pressures and densities by the 
respective values downstream. 

Figure 2.5 depicts equation (II.2) over a range of 
Macn Numbers from 1.0 to 4.0. The curve is approximated by 


the quadratic equation 


2,93 = =-2.7574 + 3.1573w - 0.2863w2 (227) 





If the high pressure is to the right the shock moves 
as Fig. 2.6 illustrates, and the governing equations are as 
listed in Table III. Again, downstream conditions are used 
to non-dimensionalize velocity, pressure, density, and the 
Riemann variable change. The development of the equations 
in Table III is given in Appendix A. A plot of equation 
(III.2) over a range of w from 1.0 to 4.0° 4s shown in yee 


2.7. The curve is approximated by the quadratic equation 


‘SrA. 2.7574 - 3.15732w + 0.2863w¢ 28 
or ( ) 
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SHOCK EQUATIONS FOR HIGH PRESSURE ON LEFT 


Definitions 
Relative Incoming Us (qn-V.) 
Mach Number: w = a Sa (iieieeesls) 
B B 
Equations 
Q Riemann 
Variable Change: 
eC OP 
aa (/+1)w ie AS a TIT y+ 
+2, Y 
7 ah) cmaw 42,1 (ie) 
i (y+L)w 
Speed of Sound Ratio: 
A 1/2 
A i oe) eee 
i, = aryl? y- 1) [1 5 W Vie 1) (II.3) 
Pressure Ratio: 
P Tier 
eee ey 2 al (IL. 4) 
P y+1 y+1 
B 
Density Ratio: 
Ze 
Pa _ _(ytl)w (II.5) 


Pp (y-1) w742 


3b 


TABLE II (CONTINUED) 


Velocity Change: 


hey Demi) 





An (y+1)w 


Entropy Change: 
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(II.6) 
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Figure 2.6 Shock Wave with High Pressure on Right 
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TABLE IIT 


SHOCK EQUATIONS FOR HIGH PRESSURE ON RIGHT 


Definitions 


Relative Incoming 
Mach Number: 








u =v. 
~~ 2s s Gage) 
An A 
Equations 
R Riemann 
Variable Change: 
= D aS A 
ey 2(1-w*) 2 B, 1 2y 2 
= Se eee LS ] + Im apy int Gay ow 
An (y+l)w  ‘y-1 A, AA ¥CyH1) Y 
2 
(C5) eo Gana) 
y+i (y+L)w 
speed of Sound Ratio: 
A 72 
2 s es ea el (ares) 
 - Temes LL pw ie q] 
Pressure Ratio: 
P al 
mee YY .¢ - yl (Gia) 
P ¥+1 v¥+1 
A 
Density Ratio: 
OQ, 2 
A. tlw (Ge 4) 
Pp (y-1)w* +2 


a2 


TABLE III (CONTINUED) 


Velocity Change: 


en Pe 2 (1=w") 


an ~  (yv+l)w 


Entropy Change: 





y-l 
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An examination of equations (II.2) and (III.2) 


reveals that 


Ey Awe 2 


This is expected since in flow to the right, the Q extended 
Riemann variable is associated with the qt+tA characteristic. 
Also, since velocity is defined as positive to the right, q 
is positive. But for flow to the left, velocity is defined 
as negative, thus q is negative. Moretti [{Ref. 3] showed 
that the downstream running Riemann variable 
(q+A characteristic) has a lower magnitude of discontinuity 
across a normal shock than the upstream running Riemann 
variable (q-A characteris-tic), and is thus preferable. [In 
flow to the left, the Riemann variable associated with 
downstream is the R Riemann variable, though the character- 
lstic associated with the R Riemann variable is q-A, in flow 


to the left gq is negative and thus 
q-A = -(iqi + 4) (2-12) 
Therefore, in fluid traveling left the R Riemann variable is 


used to find Mach Number, and in fluid traveling right the Q 


Riemann variable is used. 
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=f the shock is between two nodes with no contact 

surface within the same interval, the method to calculate 
the Mach Number, w, is as follows: 

1) The change in the appropriate extended Riemann varia- 

ble across the interval is measured; R variable for 

flow left, and Q variable for flow right. Call this 


Change, Am. 


myeeUse Am in equation (2.7) or equation (2.8) for the 
appropriate flow to calculate w. 


fe With equation (11.2) or equation (I1I.3), for flow 
right or left respectively, use w to calculate the 
analytical Riemann variable change, Ae. 

4) If this exact value of Ae is not equal to Am then 
calculate a new guess, Ag from 


[aes = OSs am = he) (2.10) 


and then repeat steps 2 to 4 until Ae equals Am. 

Thus the correct value of w is determined, and shock 
corrections from the normal shock relations are valid. V. 
is determined from equation Ciivjmeer equation (LIL. 1) 
depending on flow direction. Flow direction will be 
initially determined by the side with high pressure. High 
pressure on the left means flow will travel to the right. 
Likewise, high pressure on the right means flow will travel 
to the left. 

ee Contact Surface 

A contact surface is a boundary between regions of 
flow which are different in composition or which have 
undergone different thermodynamic histories. Thus a contact 


surface is formed when the diaphragm is burst in a shock 
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tube separating the gas initially in the driver from the gas 
initially in the driven tube [Ref. 6:pp. 23-29]. Ina wave 
rotor, a contact surface would separate the combustor gases 
from the cooler inlet air. Contact surfaces are also caused 
by two shocks of opposite families colliding or a shock over 
taking a shock of the same family [Ref. 3:pp. 15-18]. 

Figure 2.8 illustrates the changes in physical 
properties through the contact surface in the shock tube 
problem. The gas to the right has been compressed and 


heated by the shock wave, with that to the left cooled and 
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Figure 2.8 Behavior of Physical Properties Through 
a Contact Surface 


40 


expanded. The density, speed of sound, and modified entropy 
are discontinuous across the contact surface. 

Since pressure is constant across the contact 
surface, using the perfect gas equation of state and the 
first law of thermodynamics it is shown in Appendix A for 


the shock tube problem that 


(OY) (s_-s.)) 

A /A. = exp Gaya) 

Since the flow behind the contact surface to the 
tail of the rarefaction wave is isentropic, the value of Sp 
will be known. The value of Sa is the value of entropy 
behind the shock. SubEractingmeaations (174) and ~Gr.5) 
from each other and dividing by entropy Ss, Ag can be 
determined. Thus the unknown speed of sound Aq is obtained 
Meem equation (2.11). With qa equal to qa, and Sag the 
Riemann variables just behind the contact surface are calcu- 
lated from equations (1.4) and Gita )aeevOn a COncact surtace 
traveling to the left it is only necessary to interchange 
the subscripts. The velocity of the contact surface is 
easily determined since it moves at velocity gq along the q 
characteristic line associated with S. [Ref. 3:p. 16] 

4. Contact Surface/Shock Interaction 
When the shock and contact surface are within the 


same interval the calculation of the shock incoming Mach 


Number, w, must be modified. Figure 2.9 shows a plot of the 


4l 
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Figure 2.9 Modified Entropy Distribution torssncec 
and Contact Surface Traveling Right 


modified entropy change across the interval. The use of the 

Riemann values at the nodes to estimate the Riemann variable 

change across the shock would be incorrect. The correct 

Riemann variable change would be, for flow to the right, 
26723 


ASK = == (2.12) 
) B 


The correct Mach Number, w, is determined using the 
technique of Moretti modified for QAZ1D [{Ref. 3:pp. 23-24}. 
Referring to Fig. 2.9, the technique for flow to the 
right is as follows: 
1) Calculate w as if no contact surface were in the 
Interval. With values at node i+1l known, determine 
Spe Then use w and Sp in equation (11.7) to solvemaem 


Sa. This is the initial estimate of Se. Letyoe sees 


2) The Q Riemann Variable change across a contact surface 
is shown in Appendix A to be given by 
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_ yrl 
Q,1-Q51 ((S,. Sar) ( 5 )] A 


esc ca = fe Sa17Sye] x (2713 ) 
3) Define 
om oS Ades 
then 
Am = ASK + ACS (2.15) 


can be solved to obtain ASK. 


4) Using ASK as Am in the shock iteration scheme given in 
Section II.C.2 calculate a new w. 


5) With this new w and Sp use equation (II.7) to obtain a 
new Sa. 


6) Compare Sq with Sg. If they are not equal to within 
an acceptable error, calculate a new estimate for Sez 
using 


ao) (2. 16) 


Iterate until convergence is achieved. 
hits wile result in thei proper w for a condition 
with shock and contact surface interacting. 
For flow to the left, it is only necessary to inter- 


change the subscripts A and B, and use 


(y-1) 
s No ( (S Hoste ) 
- [oan ay 
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for the extended Riemann Variable change. Figure 2.10 


illustrates the left traveling condition. 
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Figure 2.10 Modified Entropy Distribution for Shee 
and Contact Surface Travelina Left 


D. BOUNDARY CONDITIONS 
1. Open Boundary 

At the open boundary, a reference pressure ratio, 
Pref = P.ofPa is specified, where P, is the pressure to be 
held at the boundary, and P, is the pressure just insidemeae 
tube. Only the case where PP, © Px 2s §eeneidercar This 
means that an outflow condition at the boundary will result 
when P. < Pa, with an expansion Wave traveling in fremm@eme 
boundary. Thus locally at the boundary, conditions are 
isentropic. 

The computational grid for the ilieft boundary is 


‘shown in Fig. 2.11. A phantom node, L, is used outside the. 
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computational mesh to enable simple enforcement of constant 
pressure and entropy at the mesh boundary node, 1. 

The approach is to determine the required variables 
at node L, then transfer these values onto node 1, and vary 
q at the boundary to meet the required boundary conditions. 
Berrst., using RBeageeac mode TL and S at node 2 in equation 


(I.3) rewritten in the form 


2 
(S Sy) (-¥(y-1))-1/¥ 
((1/P f )) 


2.18 
ref ( 


Pr ) (exp 
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Op, the density ratio at the boundary, is calculated. Using 
the perfect gas equation of state in equation (2.18) the 


temperature ratio, Tp is given by 


2 
(y (ley) (S a ; 


eee) 2 ee 
T, = (0, ) (exp ( ) 
and can be calculated once pp is known. With an initial 


velocity at the boundary, qp, assumed the Riemann variables 


R and Q are determined for node L using 


R = q,- "Tp s (22a) 

and 
Oe Vicia a VTP S | (2.21) 
Subtracting equation (e200) from ( 2ean2 aa and 


dividing by twice S yields A at node L, where A = vyRcT. 
The values of R, Q, A, Gp, and S are now Known at nodeuie 
These values are assigned to node 1, and the QAZ1D algorithm 
is applied to the boundary node as if it were an interior 
node. The result is an estimate of the new values for R, Q, 
and S at node 1. These are used to obtain a new estimate of 
the pressure ratio at node 1, call it Ppy. Ppn is compared 
with Pros, and the new value of S with the old value Oreage 


node 1. A new qp is calculated by solving equations (2.20) 
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and (2.21) simultaneously. The process is repeated until 
entropy and pressure remain constant at the boundary. For 


the right-hand boundary, Fig. 2.12 shows the computational 


a 
iy 


right grid boundary 
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Figure 2.12 Right Boundary Computational Grid 


grid with the phantom node, r, to the right of the mesh 
boundary node, n. The procedure for the right boundary is 
the same as for the left boundary. 

When a shock travels across the open boundary, the 
boundary is first corrected using the method described in 
Beecion I1.C.2. Subsequently, the pressure is returned to 


Give the constant value specified for Pyar, using the above 
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procedure. The result will be an expansion wave traveling 
inward. 

For supersonic flow, all of the variables are deter- 
mined from values at the interior nodes [Ref. l:p. 3]. 

2. Closed Boundary 

The solid boundary is imposed by setting the 
velocity at the boundary, gq, equal to zero. The same 
computational grid is used as that for an open boundary, but 
a different technique is used in assigning values to the 
phantom node. Referring to Fig. 2.11, the velocity at node 


2, G9, is known. By setting 


dy, = - q (2.925 


the wave impacting the boundary will be met by a wave of 
equal velocity but opposite in direction and will result in 
node 1 appearing as a solid boundary. Further, Wea 


facilitate the QAZ1D method, set 


Rr = ae (2.23) 
Q, = |Rol (2 s285) 
Sr = So (2425 
A = .2°A (2.26) 
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This will enable the computation of the new boundary node 1 


values as an interior point with 
q = 0 


The right boundary is handled in the same way. 
Figure 2.13 illustrates the sequence of events when 


a shock reflects off a solid boundary on the right [Ref. 7: 





=O 
as 
sod solid 
boundary Be boundary 
1 i) |e 
A) Shock Approaches Solid B) Shock Reflects off 
Boundary Solid. Boundary 


Figure 2.13 Shock Reflecting at Solid Boundary 


Pee 172-195). Since the velocity behind the reflected 
shock, gp, is zero, and the velocity in front of the shock, 
Ga, 1S known, the velocity gradient, Aq, over the shock is, 


in non-dimensionalized form, given by 
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Aq = (qg - Ga)/Aqg (2°. 259 


Rearranging equation (III.6), and using Y= 1.4 it 


can be shown that 
7 
w = 1 +y0.36\(Aq) > = Oem (2.289 


This is an exact analytical value for w based on stationary 
normal shock relations. The speed of sound ratio, Ap/Aq, is 


calculated from equation (III.3) and 


Ap = Ag(AB/Ag) (2.208 


Entropy behind the shock, Sp, is determined from) equauua 
(III.7) since Sa, 1s known. Then the extended Riemann 
variables, R and Q, behind the reflected shock are 
calculated from equations (1.4) and (1.5) respectively. 

For a shock reflected from the left boundary, the 
subscripts A and B are interchanged and using equation 


(II.6) the Mach Number is given by 
w = V1 + 0.36(Aq)2 + 0.6iq (2 .2ap 


The equations from Table II are used instead of 
those from Table III (used in the case of reflection from 
the right boundary), since the shock is now traveling from 
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Ziz=. FORTRAN PROGRAM E1DV2 


A. GENERAL DESCRIPTION 

Program E1DV2 is a Fortran computer program which 
calculates l1-dimensional unsteady flow based on the QAZ1D 
formulation and method of solution of the Euler equations of 
motion. The program has provisions for tracking in time the 
location of shock waves, contact surfaces and head and tail 
of the expansion waves. Its intended application to wave- 
rotor design and analysis explains the present definition of 
node points and boundary conditions. tnawes CULrent form, 
the code describes flow in a tube initiated by a diaphragm 
bursting. The location of the diaphragm, whether the flow 
is to the left or the right, and whether the ends of the 
tube are closed or open, can be varied. 

The E1DV2 program is written in FORTRAN 77, and was 
developed using the IBM 3033 System 370 computer. The use 
of DISSPLA subroutines in plotting results requires 
compiling and executing using VS FORTRAN procedures. 
Appendix B describes the procedures for running the E1DV2 
program on the Naval Postgraduate School IBM computer 
system. A listing of the code, which contains its own 


description in comment statements, is given in Appendix D. 
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The code incorporates: 
1) first order accuracy 


2) double precision numerics, except for single precision 
in graphical output 


3) linear interpolation and extrapolation algorithms 

4) quadratic polynomial approximations for curve fitting 
5) non-dimensional variable input and output 

6) structured subroutine format. 

Table IV lists all the subroutines called from the main 
program and brief descriptions of their purposes. The 
extensive use of subroutines was intended to allow modifica- 
tions and extensions ole the code without major 


restructuring. 


B. THE MAIN PROGRAM 

The main program flowchart is illustrated in Appendix C, 
Eig, “Cale: The conventions and a comprehensive list of 
variables are listed in the beginning of the program. The 
number of grid nodes is set by the user, and must be an odd 
number. Arrays are dimensioned to the number of nodes in 
the main program prior to compiling and executing. 

Initial values for temperature, pressure, and density 
ratios at the diaphragm and boundaries are entered by 
editing. Also, the initial velocity distribution on each 
side of the diaphragm is set, and as well as the ratio of 
specific heats. The side with the high pressure is identi- 


fied. Boundaries are set either closed or open. In the 
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TIME 


TRAK 


SWEEP 


COND1 


COND2 


COND3 


COND4 


CONDS5 


COND6 


COND7 


COND7N 


TABLE IV 


SUBROUTINES 


Calculates minimum time step 


Tracks discontinuities, calculates shock mach 
number, w 


Advances node values to next time step 


Calculates interpolated values of Riemann 
variables, 3q/ds, and dA/ds when no shock 
or contact surface within an interval 


Calculates interpolated values of Riemann 
Variable, %/3s, and dA/ds when discontinuity 
within interval to right of node or flow is 
supersonic to the right 


Calculates interpolated values of Riemann 
varlable, “/3s, and dA/ds when discontinuity 
within interval to left of node or flow is 
supersonme to the left 


Called when node is jumped by discontinuity 
to load node information into array for use 
in subroutine CORRCT 


Calculates interpolated values of Riemann 
variable, j9q/3s, and dA/ds when discontinuity 
on either side of node and flow is 

supersonic 


Would contain subroutine to advance node 
when contact surface is in front of shock 
after intersecting and traveling in same 
direction. Not currently in use 


Would contain subroutine to solve precisely 
when and where shock and contact surface 
would intersect within an interval. Not 
currently in use 


Would ccntain subroutine to solve precisely 
when and where shock and contact surface would 
intersect if either were jumping a node simul- 
taneously. Not currently in use 
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COND7S 


COND8 


COnRew 


DELTAX 


SKJOUMP 


CSJUMP 


EXTRAP 


TNR 


DBURSE 


BONDRY 


SRE LCL 


“BBD Y 


BORDER 


PLOT 


EXACT 


bier 


TABLE IV (CONTINUSD) 


Would contain the subroutine to solve the 
Riemann problem of shock and contact surface 
intersecting. Not currently in use 

Would contain subroutine to advance node after 
shock and contact surface have crossed and are 
moving apart. Not currently in use 


Advances nodes jumped by discontinuity to next 
time step 


Determines relative location of discontinuity 
within an interval 


Calculates variable change across a normal 
stationary shock 


Calculates variable change across a contact 
surface 


Extrapolates entered values 
Interpolates entered values 


Calculates initial Riemann values at node 
where diaphragm is located 


Calculates values to update left and right 
boundary 


Calculates new shock parameters when shock 
reflects off solid boundary 


Alternate interpolating scheme near boundary 
Sets graphical borders 

Plots gq, S, FF, andae 

Plots exact 0 versus calculated 0 


Lists calculated values in files 8, 9, 10 
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open case, either constant pressure, or an adjustable 
pressure is specified. The latter case in effect extends 
the tube, and allows discontinuities to disappear out of the 
tube without creating expansion or compresSion waves. Exact 
values for velocities of expansion wave, shock, and contact 
surface, and density ratio are specified by the user in the 
input section. 

Output is controlled by setting the variable GRAPHS to 
either 0, 1, or 2. A zero creates three files which contain 
data calculated at the time the call to subroutine LIST is 
made. A value of 1 causes plots of pressure, density, 
velocity, and entropy distributions to be created. When 
GRAPHS equals 2 a plot of the exact density distribution 
along with the calculated density distribution is made. 
Carls to output routines are determined by the parameter, 
SKIP, the number of time steps between calls, which is 
specified by the user. 

Program termination can occur for several eAsene 
First, the program terminates when a maximum number of time 
steps, JSTOP, is reached. JSTOP is specified by the user. 
second, the program terminates if, during execution, any of 
the following conditions are met: 

1) The shock intersects with the contact surface 


2) The pressure at any open boundary is higher than the 
pressure at the first node inside the boundary 


3) The shock exits an open boundary. 
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In these cases a message is displayed on the terminal 


describing the reason for termination. 


C. THE SUBROUTINE PROGRAMS 
1. The "DBURST" Routine 

Figure C.2 in Appendix C shows the flowchart for 
Meno Ls The assumption that a shock and contact surface 
are formed immediately causes oscillating numerical solu- 
tions that dampen within five to ten time steps. By 
mathematically "bursting" the diaphragm prior to time zero, 
a solution at the node where the diaphragm is located can be 
determined at time zero. The result is a more accurate 
representation of the wave structure and a dampening of a 
reduced transient solution within three time steps. 
"DBURST" calculates the Riemann variables that would exist 
behind a shock and contact surface that have moved an 
infinitesimal amount after the diaphragm burst, and assigns 
them to the node where the diaphragm was situated. 

2. Lhe OIE! Routine 

The "TIME" subroutine was not changed from that 
developed by Salacka [{Ref. 2:p. 35]. The purpose is to 
compute the maximum allowable time step but ensuring that 
all characteristic trajectories over the time step are 
within one interval between nodes. The time step is 


determined, from 
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DELT = H/ABS[QH] (3.1) 


at every node, and the minimum time step is used. 
3. The "TRAK" Routine 

This routine computes the new locations of the 
shock, contact surface, and head and tail of the expansion 
wave after the time step, DELT, calculated in "TIME." £In 
memeendix C, Fig. C.3 shows the flowchart for the "TRARK" 
subroutine. The shock velocity, Ve, incoming Mach Number, 
w, and pressure, density, and sonic velocity ratios PR, DR, 
and AR respectively are determined. The equations for flow 
left and right are the same, except for DQ, the velocity 
gradient, which is different only by sign. Thus by proper 
bookkeeping the same equations are coded once and used in 
either direction. The initial direction of the shock and 
contact surface are set from code that interprets which side 
has the high pressure at time zero. The contact surface 
speed is determined after the shock speed is established. 
The situation at time zero is handled as in "DBURST" to 
ensure that the shock/contact surface interaction is 
properly modeled. Since there is a certain transient 
solution that decays after three time steps, tracking of the 
expansion wave is not initiated until then. 

The head of the expansion wave travels at a velocity 
corresponding to qtA, and the tail with speed q-A. In the 


case of flow traveling to the left, this is reversed. 


7) ° 


4. The “SWEEP! Reueane 

The "SWEEP" subroutine solves equation (2.5) and 
updates the values of variables at the nodes, including the 
boundary nodes. Appendix C, Fig. C.4 auete the "SWEEP" 
Eiewenawt-. Only those nodes which have been crossed by a 
discontinuity during the time step are updated by the 
"CORRCI" subroutine. Section II.B describes the theory 
coded into "SWEEP", The routine begins at the left boundary 
node and progresses to the right boundary node, one node at 
a time. Figure 3.1 illustrates the overall algorithm 


applying the QAZ1D method. At the boundary nodes, the 





Figure 3.1 Overall Algorithm for QAZ1D Method 
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"BONDRY" subroutine is called and returns the updated values 
at those nodes. In between, the correct and most effective 
algorithm to update the node is determined. These 
algorithms are coded into eight different condition 
subroutines. These routines calculate the interpolated 
Riemann variables associated with the three characteristics 
(q+A, q, and qA), plus the spatial derivatives at dq/3s and 
aA/aSs. This allows "SWEEP" to calculate Aw and 2z. The 
integral of z is then determined, and the update for the 
variables of that node in dw are computed. Ww is stored 
until the entire mesh has been swept. Then the variable 
arrays (W) are updated to the next time interval. 

To choose the proper condition routine, the 
following information about conditions in the interval on 
Sener Side of the node are expressed in two 3 digit integer 
variables, SHOCK and CNTACT. The value of these variables 
provide a code for the following information: 

SHOCK--The existence of a shock within an interval, and, 
if so, the direction of travel, and if the shock 
will cross the node in front of it within the 
current time step. | 

CNTACT-The existence of a contact surface within an 
interval, and, if so, the direction of travel, and 
if the contact surface will cross the node in front 
of it within the current time step. 

Figure 3.2 illustrates the regions around a node. Table V 
lists the values that SHOCK and CNTACT may have, and their 


meaning. Figures 3.3, 3.4, and 3.5 illustrate examples of 


what different SHOCK and CNTACT values mean. At each node 
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Figure 3.2 Interval and Node Description used in "SWEEP" 
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Figure 3.3 Schematic of SHOCK = 331 and CNTACT = 232 
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TABLE V 


VALUES OF PARAMETERS SHOCK AND CNTACT 
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SHOCK 
Located Direction Crosses node in 
Value in interval on of travel path within time step 
100 No Shock N/A N/A 
2 1) Left Right Yes 
322 Left Right No 
eer 1 Left Left Yes 
332 Left Left No 
B21 Right Right Yes 
222 Right Right No 
23 1 Right Left Yes 
232 Right Left No 
CNTACT 
100 No contact surface N/A N/A 
B21 Left Right Yes 
Bie 2 Left Right No 
Bo 1 Left Left Yes 
332 Left Left No 
Ze |. Right Reght Yes 
ee 2 Rect Right No 
23) SULe gue Left Yes 
eo Right Left No 
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Figure 3.4 Schematic of SHOCK = 100 and CONTACT = 321 
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Figure 3.5 Schematic of SHOCK = 221 and CNTACT = 222 


the code examines the shock and contact surface condition, 
along with relative location of the shock to a contact 
surface if both exist near the node, if the flow is super- 
sonic or subsonic at the node, and the direction the flow is 
traveling. The appropriate algorithm is determined and 
called in the form of a condition subroutine. 

The "SWEEP" routine is designed to handle all 
possible shock and contact surface combinations. Because 
the code does not solve the situation where the shock and 
contact surface intersect, as illustrated in Fig. 3.6, then 
any combinations after the shock and contact surface have 
crossed call condition routines that currently contain only 
messages. The code was intentionally prepared to be easily 
extended to treat the intersection of the shock and the 


contact surface, and beyond. 


shock and 
contact surface contact surface snock 
at t —— >> intersectat t+dt qW att 


q 





Figure 3.6 Schematic of SHOCK = 231 and CNTACT = 322 
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5. The "Condition" Subroutines 
The ten subroutines that are called by "SWEEP" to 
calculate the interpolated values of the extended Riemann 
variables, R and Q, plus S, dq/ds, and dA/d3s are CONDI, 
COND2, COND3, COND4, CONDS5, COND6, COND7, COND7S, COND7N, 
and CONDS8. The procedure used in COND1, COND2, COND3, and 
COND5 is that of Salacka [Ref. 2:pp. 36-37] modified to 
account for contact surfaces. Appendix C, Fig. C.5 ism 
general flowchart for these four routines. Essentially, an 
initial estimate is made of the characteristic slopes, }i, by 
enforcing the principle of domain of dependence. The 
assumption is made that the slopes are linear, and that the 
Characteristics pass through point B, as defined in Fig. 
Sy Then q and A are computed at point A, which in turn 
yields a second estimate of the characteristic slope, }. 
The two slopes for each characteristic are compared. 
If they are not equal to within an acceptable error, the new 
estimate oof \} 1S used to repeat the process’ until 
convergence is achieved. All three characteristics are 


handled simultaneously. The value of 
Ns = owe (32 


is then determined for each characteristic qt+A, q, and q-A. 
This allows linear interpolation of Q, R, and S at point A. 
The spatial derivatives, 30/3s and 9A/0S are computed from 
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Figure 3.7 Computational Method for COND1 Routine 


the values of q and A used in the estimate of the initial 
characteristic slopes. 

The "COND1" subroutine is used when there are no 
contact surfaces or shocks within the interval on either 
Side of the node, and flow is subsonic. The qtA character- 
istic uses forward differencing, while the q-A characteristic 
uses a backward differencing scheme to keep the initial 
domain of dependence of the numerical scheme outside the 
physical domain of dependence. 

The "COND2" subroutine is a backward differencing 
algorithm used in situations when one or more discontinui- 
ties are in the right interval, or if flow is supersonic to 


the right. Fig. 3.8 illustrates the computational method 
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Figure 3.8 Computational Method for COND2 Routine 


for “COND2: The characteristic slopes of q+A and q are 
handled as in "COND1", but the q-A characteristic slope is 
determined by a backward rather than a forward scheme. To 
interpolate between node i and i+l would be incorrect 
because of the discontinuity in the values between them. 
The value of parameters in.the right interval are interpo- 
lated from values in the left interval by assuming the 
derivatives do not change between the intervals up to the 
Giscontinuities (i.e., a shock and/or a contact surface). 
The "COND3" -subroutine is a forward differencing 
algorithm used in situations when one or more discontinui- 
ties are in the left interval, or if flow is supersonic 


traveling to the left. Figure 3.9 shows the computational 
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Figure 3.9 Computational Method for COND3 Routine 


method for "COND3". The method is similar to that in COND2 
but the qH characteristic slope has to be interpolated from 
values of node i instead of node i-1 using a forward 
difference scheme. The characteristic slope for q-A and gq 
are estimated by forward difference schemes using values of 
node i+1 and i respectively. 

For conditions, such as that illustrated in Fig. 
3.10, where a discontinuity is in the interval opposite the 
direction of supersonic flow then COND5 subroutine is 
called. A mix of COND2 and COND3 is used. For flow to the 
right the method of COND2 is implemented, while flow to the 


left is the same as for COND3. 
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Figure 3.10 Example Condition fer "CONDS”™ Rover 


Subroutine "COND4" is called whenever a node is to 
be crossed by a discontinuity from either direction. The 
flowchart for "COND4"is shown in Appendix C, Fig. C.6. ‘Two 
integer vector arrays, LNODE and RNODE, are used to store 
the following information on the node being crossed. 

A) the node number, I 

B) the value of SHOCK 

C) the value of CNTACT 

D) the current value of the integer time step, J 
The information is used in "CORRCT" to update the node to 
the next time level. In the current code the maximum number 
of nodes that can be crossed during any time step is two. 
Figure 3.11 shows an example of this situation, and defines 
the assignment of values. Because the nodes are swept from 
left to right, the left-most node is assigned to LNODE, 
while RNODE applies to any node crossed to the right of the 
LNODE node. The values of AW are set to zero, so 
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Figure 3.11 Example of "COND4" Routine Situation 


that these nodes are skipped when updating occurs in 
SWEEP". The main program interprets the value of J in 
LNODE and RNODE to determine if zero, one, or two nodes need 
m@emebe corrected in "CORRCT". If zero nodes are crossed 
during a time step then "CORRCT" is not called. 

The "COND6" and "COND8" subroutines are designed for 
future extensions of the entire code. Until the Riemann 
problem illustrated in Fig. 3.6 1s coded, then any situation 
after the shock and contact surface intersect cannot be 
computed. However, "SWEEP" is already coded to call "COND6" 
fer ene situation illustrated in Fig. 3.12, or its mirror 
image. "“COND8" would be called for all the other situations 
after the shock and contact surface interaction has taken 
place. Currently, both subroutines output messages on the 
screen describing the situation, and set HALT equal to 1 to 


Germainate the program. _ 
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Figure 3.12 Example of "COND6" Situation 


The "COND7", "COND7N", and "COND7S" subroutinesmana. 
all related to each other. They are intended to contain 
code that would facilitate the shock and contact surface 
intersection and solve the resulting Riemann problem. The 
"COND7" routine is called when the contact surface and shock 
moving in opposite directions would cross in a time step 
within an interval. The routine would calculate when and 
where within the time and space domain the intersection 
would occur based on the known velocity of each discontinu- 
ity and their current locations. This new time could then 
be used to rerun "SWEEP" with the shock and contact surface 
exactly on top of each other at the end of the new time 
step. The "COND7N" routine is for the same situation, but 
when a node is crossed by either discontinuity during the 
same time step. The same procedure is done, with the 
requirement to ensure the node crossed 1S properly updated 
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by “CORRCT" re new time step. The "COND7S" routine 
is called when a shock and contact surface are located at 
the same point at a time other than zero. The subroutine 
would contain the code to solve the Riemann problem set up 
by "COND7N" or "COND7" routines. 

6. The "CORRCT" Routine 


Appendix cC, Fig. cC.7 shows the flowchart for 


eeRRCT", The "CORRCT" routine corrects nodes that have 
been crossed by a-discontinuity (1i.e., a shock or a contact 
surface) or are straddled by both. The routine takes the 


information from LNODE and RNODE arrays to determine which 
node or nodes are to be updated, and if there is any shock 
and contact surface interaction at the nodes. Then, using 
subroutines "DERGAX "SKJUMP", TeCou UMP, BS Neen 
"INTERP", and "BBDRY" the concepts from Section II are 
imposed to calculate the jump in the parameters R, Q, S,A 
and q at the node in question to the new time level. 
7. The "BONDRY" Routine 
The "BONDRY" routine computes the new values of the 
parameters Q, R, and S at the boundary nodes for time step 
Ges The concepts of Section II for an open or closed 
boundary are coded as shown in the flowchart of "BONDRY" 
given in Appendix Cc, Fig. C.8. The routine uses the 
weenor", "“COND2", or "COND3" routines to calculate the 
correct w, and then solves for w the same as in "SWEEP". 


If a shock crosses an open boundary node during a time step, 


efit 


then "SKJUMP" routine is used to caiculate the values behind 
the shock at the node. 
8. The "SRFLCT" Routine 

The "SRFLCT" subroutine is called to calculate the 
values of Q, R, S, q and A at the boundary node when a shock 
reflects from a solid boundary. Appendix C, Fig. C.9 shows 
the flowchart of "SRFLCT" which codes the analysis in 
Section II.D.2 for closed boundary shock reflection. The 
time for the shock to reach the solid boundary, oty, is 
computed. Then the excess time in the time step, ot, is 


given by 


After the new shock velocity, V.s,is calculated, then the Hew 
location of the reflected shock is known from multiplying V. 
by dtey and adding or subtracting from the boundary 
location. 
9. The "DELTAX" Routine 
The "DELTAX" subroutine is used in "“CORRCT" and 
"BONDRY" to calculate the location within an interval of a 
discontinuity in terms of x. Figure 3.13 defines the 
various displacements which are calculated. 
102 The "SKIUMP! Roucemne 
The "SKJUMP" subroutine computes the conditions 


downstream of a stationary normal shock as described in 


ee 


discontinuity 


eft 
undary x? 





Figure 3.13 Discontinuity Location within an Interval 


Section 2 ig) Gs B ge timer “CORRCT"™, RAK, and "BONDRY" 
subroutines call "SKJUMP" when required, to obtain the speed 
of sound (A), entropy (S), and velocity (q) behind the 
shock. | 
11. The "CSJUMP" Routine 
The "CSJUMP"" subroutine computes the velocity and 
speed of sound change across a contract surface as described 
mmeeoection I1I.D.3. The "CORRCT" and "TRAK" routines call 
eoUUMP". 
12. Numerical Support Routines 
The "EXTRAP", "“INTERP", and "BBDRY" routines are 
oeoy "CORRCT"’, "“BONDRY", "TRAK", and "“DBURST" to extrapo- 
late or interpolate data to the face of a discontinuity or 
Beant. in particular, "BBDRY" is an interpolation routine 


used within an interval near a boundary. 
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13. The Output Routines 


The four output subroutines used are those developed 
by Salacka [Refs 2p. 39-40)- The "BORDER" and "PLOT" 
routines output plots of pressure, velocity, density, and 
modified entropy distributions versus a non-dimension- 
alized tube length at preset time intervals. "EXAGi 
creates a plot of the exact density distribution at selected 
points and compares it with the computed density 
distribution for a selected time step. The "LIST" routine 
outputs to files 8, 9, and 10 tabular listing of computed 


data. The files are sent to the user's permanent storage 


disk. 

A value of GRAPHS equal to one causes "BORDER" to 
be called in the main program. This defines plot axis, 
labels, and headings. "PLOT" is then called once at time 


zero, and at every time step set by SKIP. 

If GRAPHS is set equal to two, then "EXACT" is 
called every SKIP time steps to plot six exact density 
values and the computed density distribution. The six exact 
points are: 

A) The two boundary points 
B) The head and tail of the rarefaction wave 
C) A point just behind the contact surface 


D) A point gust behindsehe= shock 
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The location of the exact values is based on elapsed time 
and known exact values for wave velocities entered with the 
initial conditions for the problem. 

"LIST" routine is called every SKIP time steps when 
GRAPHS is equal to zero. File 9 contains a tabular listing 
of the Riemann variables, modified entropy, pressure, 
density and velocity distributions, elapsed time, discon- 
tinuity velocity and location. File 10 is a tabular listing 
of the location of the shock, contact surface, and the head 
and tail of the expansion wave computed at every SKIP time 
steps. File 8 contains the exact values of the location of 
the shock, contact surface, and the head and tail of the 


expansion wave. 
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IV. RESULTS 


Four test cases were run using the E1DV2 code. The 
purpose was to verify the ability of the code to determine 
the unsteady flow process and correctly simulate varying 
boundary conditions and flow directions in the Riemann shock 
tube problem. The shock tube is illustrated in Fig. 4.1 
with the high pressure on the left. The tube is divided 
into two sections by a diaphragm. When the diaphragm is 
burst, the pressure equalizes through a shock wave traveling 


into the expansion chamber, and an expansion or rarefaction 


diaphragm 


r expansion chamber 
high 


compression chamber} 


expansion wave contact 
; shock 





time =t 


Figure 4.1 Shock Tube at t = 0 andt=t 


wave traveling into the compression chamber. A contact 
discontinuity is behind the shock wave, traveling at the 


particle velocity [Ref. 6:p. 30). 


meee test CASE 1 

Test Case 1 was designed to demonstrate shock reflection 
and expansion wave reflection at solid boundaries. Test 
Case 1 had the following initial and boundary conditions: 


1) Pressure and density ratios equal to five, with high 
pressure on the left 


2) Temperature ratio equal to unity 

3) Both boundaries were closed 

4) Diaphragm was located at x = 0.5 

5) Computational mesh had 101 nodes 

6) Maximum time step, JSTOP, = 109, with SKIP = 18. 
Plots of the results obtained for the pressure, density, 
velocity, and modified entropy distributions are shown in 
Meee 4.2 and Fig. 4.3. Fig. 4.2 shows the computation up to 
time step 55, while Fig. 4.3 takes the computation from time 
feeoe>6o to termination. The output of "EXACT", a plot of 
the exact density compared with the computed density 
distribution is shown in Fig. 4.4. Fig. 4.5 illustrates the 
spatial location versus time for exact and computed shock, 
contact surface, and head and tail of the expansion wave. 
Data for Fig. 4.5 were taken from file 10 and file 8. Fig. 
4.6 shows plots of pressure, density, modified entropy, and 


velocity distributions for the same initial and boundary 
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SHOCK TUBE RESULTS 
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Figure 4.2 Test Case 1 (J = 1 to J = 55) 
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Figure 4.6 Test Case 1, High Pressure on Right 


82 


conditions but with high pressure on the right. Figure 4.6 
includes data corresponding to those in both Fig. 4.2 and 


Pig. 4.3 for the high pressure initially on the left. 


B. TEST CASE 2 

Test Case 2 was to demonstrate an open boundary, 
constant pressure, expansion wave interaction. The test 
case was run with the following initial and boundary 
conditions: 


1) Pressure and density ratios equal to 5.0, with high 
pressure on the left 


2) Temperature ratio equal to unity 


3) The left boundary was open, with a constant pressure 
ratio across the boundary of 4/5 


4) The right boundary was open, with an adjustable pres- 
sure ratio that in effect extended the length of the 
tube 

5) The diaphragm was located at x = 0.74 

6) Computational mesh had 51 nodes 

7) Maximum time step, JSTOP, = 37, SKIP = 9. 

Plots of the results for the pressure, density, velocity, 
and modified entropy distributions are shown in Fig. 4.7. 
The test case was rerun with high pressure on the right. 
Figure 4.8 shows the results for pressure, density, and 
temperature ratios the same. The diaphragm was then at x = 


0.26. Boundary conditions were reversed. The same computa- 


tional mesh, and output control were used. 
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Figure 4.7 Test Case 2, High Pressure on Left 
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Figure 4.8 Test Case 2, High Pressure on Right 


85 


CC. WTEST CASE 3 

Test Case 3 was designed to demonstrate a shock exiting 
an open boundary with constant pressure maintained at the 
boundary itself. The following initial and boundary 
conditions were set: 


1) Pressure and density ratios equal to 5.5, with high 
pressure on the left 


2) Temperature ratio equal to unity 
3) Left boundary was closed 


4) Right boundary was open with a constant pressure ratio 
across the boundary of unity 


5) Diaphragm was located at x = 0.5 

6) Computational mesh had 101 nodes 

7) Maximum time step, JSTOP, = 73, with SKIP = 18. 
Figure 4.9 shows plots of the results obtained for pressure, 
density, modified entropy, and velocity distributions. 
Similarly, Fig. 4.10 shows results obtained by putting the 
high pressure on the right, reversing the boundary condi- 


tions, and holding everything else the same. 


D. TEST CASE 4 

Test Case 4 was designed to demonstrate a lower pressure 
ratio, with shock wave reflection. The initial and boundary 
conditions were set as follows: 


1) Pressure and density ratios equal to 3.2, with high 
pressure on the left 


2) Temperature ratio equal to unity 


3) Both boundaries were closed 
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Figure 4.9 Test Case 3, High Pressure on Left 
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Figure 4.10 Test Case 3, High.Pressure on Right 


- Bis 


4) Diaphragm was located at x = 0.5 

5) Computational mesh had 901 nodes 

6) Maximum time step, JSTOP, = 601, with SKIP = 150. 
Plots of the results obtained for the pressure, density, 
modified entropy, and velocity distributions are given in 
Fig. 4.11. Figure 4.12 shows the results obtained with the 
high pressure to the right, SKIP = 200, and holding all 


other parameters constant. 
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Figure 4.12 Test Case 4, High Pressure on Right 
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V. DISCUSSION 


A. RESULTS OF Tisai] cases 

Overall, in the four specific test cases which were run, 
the expected qualitative flow behavior was produced by the 
code. Quantitatively, either experimental data or further 
exact solutions are needed to fully validate the computa- 
tions. The results of each test case will be discussed 
separately. 

Test Case 1 results in Fig. 4.2 show the formation of a 
well-defined shock wave traveling to the right, with the 
contact surface crisply defined following along behind. The 
entropy drops sharply across the shock and remains constant 
to the contact surface and then jumps to a steady, constant 
value to the left boundary as expected. Pressure and 
velocity remain constant through the contact surface. 
Velocity is positive since flow is to the right. The 
expansion wave is smeared from head to tail over the correct 
range. Figure 4.3 is a continuation of the flow problem 
with a sequential display of the results. The entropy 
continues to drop as the reflected shock travels back toward 
the contact surface. Notice that the velocity behind the 
shock is zero. At the left boundary as the head of the 
expansion wave reflects, pressure and density continues to 


drop while velccity 1S zero at the boundary. The density 


SZ 


distribution shows that the contact surface and shock are 
about to cross. The program terminated and issued a message 
before the next output call was made, when the _ code 
determined that the contact surface and shock would cross. 
The exact density distribution compared with the computed 
density distribution in Fig. 4.4 clearly shows the shock and 
contact surface are modeled very accurately with a medium 
mesh. The expansion process, consisting of infinitesimal 
changes propagating along the characteristics is fairly well 
modeled. Figure 4.5 shows that tracking of the contact 
surface and shock wave, in comparison with exact locations at 
set times, is excellent. Using q+A and q-A based on condi- 
tions behind the burst diaphragm as estimates of the 
velocities of the head and tail of the expansion wave 
respectively results in positions of the head and tail 
Slightly different from those of an exact Riemann solver 
code [Ref. 8]. The exact code is based on a method given in 
Peet. O9°pp. 181-191). However, an examination of the 
tabulated output showed that the solution method predicted 
epeniges in the gradients to occur at similar locations to 
those computed using the qt+tA and q-A estimates. Consequent- 
ly, while there may be some inaccuracy in the numerical 
solution, the method of tracking the Mesa and tail of the 
expansion wave is acceptable. 

The ability to reproduce the same conditions but with 


flow to the left is demonstrated clearly in Fig. 4.6. The 


oe 


entropy, density, velocity, and pressure distributions are 
mirror-images of those for the case of high pressure on the 
left. 

Test Case 2 results in Fig. 4.7 show that at time zero 
with a pressure ratio at the boundary of 4, which is lower 
than the preset value inside the tube of 5, an expansion 
wave traveling inwards is generated. Outflow is seen in the 
negative velocity distribution developing near the left 
boundary, where negative velocity implies flow traveling to 
the left. The contact surface and shock are clearly formed 
to the right of the off-center diaphragm, and travel to the 
Pui, The expansion wave generated at the diaphragm 
travels to the left and intersects with the expansion wave 
generated at the left boundary. The pressure distribution 
shows the propagation of these waves as the pressure drops 
behind the head of each wave. With conditions reversed, so 
that the high pressure is to the right of the diaphragm, the 
results in Fig. 4.8 are a mirror-image with velocity 
negative for flow to the left, and positive for the outflow 
COuvene: cage. 

Test Case 3 results in Fig. 4.9 demonstrate that the 
shock wave meeting an open boundary where the pressure is 
heid constant, 15 correctly, Sonanted to exit the tube. The 
density distribution shows the jump increase across the 
shock, then another jump increase across the contact 


surface. Then density plunges down as an expansion wave 
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travels back into the tube when the shock exits. The 
pressure ratio at the open right boundary was held at unity, 
as seen in the pressure distribution which, at the boundary 
behaves similarly to the density distribution. The pressure 
at the boundary would be required to increase as the shock 
passed by, however the enforcement of the constant pressure 
locally causes an expansion wave to form traveling to the 
left. The entropy remains constant, at the value behind the 
shock, from the boundary back to the contact surface. There 
is a jump in entropy across the contact surface. The 
velocity distribution shows that as the shock travels to the 
Beene, the velocity jumps up. When the expansion wave 
generated at the right boundary travels inward the velocity 
continues to increase in magnitude while flowing to the 
right. The expansion wave generated at the diaphragm can be 
seen traveling to the left in the plots of pressure, 
density, and velocity. Reversing the situation and 
Semputing conditions for flow to the left, in Fig. 4.10, 
results in a mirror-image in the pressure, density, and 
entropy distributions. The velocity becomes negative since 
flow is to the left. 

Test Case 4 is a repeat of the first test case but with 
a pressure ratio of 3.2 and a fine computational mesh. In 
Fig. 4.11 the density steps up across the shock and contact 
surface as they travel to the right. When the shock 


reflects the density jumps again. The velocity distribution 
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shows that the velocity drops to zero behind the reflected 
shock, after it had jumped up across the shock when it was 
traveling to the right. The instability seen in the 
velocity and pressure distributions near the midpoint 
occurred during the first 100 @timessccepa. A numerical 
instability appears to generate at the diaphragm at time 
zero for pressure ratios less than about 3.0 which can be 
severe. For the conditions in this test case the transient 
instability damped out after the first 100 time steps. The 
shock, contact surface, and expansion wave are nevertheless 
seen to be accurately modeled. Figure 4.12 demonstrates 
that reversed conditions result in mirror-images of the 


computed conditions. 


B. CURRENT LIMITATIONS OF THE CCDE E1DV2 

The current limitations of the E1DV2 code for modeling 
quasi-one-dimensional inviscid flow are in two categories. 
First, there are numerical limitations in obtaining 
solutions with the present code for different initial and 
boundary conditions. Second, the present coding limits what 
flow situations can be treated. 

The numerical instability which occurs at low pressure 
ratios during the first few time steps can be severe. The 
assumption in the current code is that the shock forms 
immediately, which at high pressure ratios does not pose any 
problems. The mathematical modeling of the bursting 


diaphragm in subroutine "DBURST" adequately reduces the 
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teansient instability for high pressure ratios, but not for 
low pressure ratios. 

A second rather different numerical limitation was 
identified in numerically feteceing the location of the head 
and tail of the expansion process. Because the expansion 
process is not a sharply defined front, fixing the precise 
locations of the head and tail waves numerically at a given 
time is difficult. However, the method currently used does 
accurately track the computed propagation of the wave along 
the characteristics. The characteristics are modeled 
currently as straight lines. A higher order curve would 
describe them more accurately. 


The current code is limited to tracking two discontinui- 


ties and the expansion wave (i.e., a shock and a contact 
surface). Thus any situation which would generate another 
discontinuity can not be computed. The code can determine 


when this will occur, and will issue a message explaining 
the situation before terminating. Thus the shock colliding 
with a contact surface will currently terminate the program. 
A boundary pressure that is higher than the pressure inside 
the tube would also create a compression wave or possibly a 
shock wave. This condition will also terminate the program. 
Simulation of the process of a shock forming over an 
extended distance and time as compression waves pile on top 
of each other and: strengthen, would also require additional 


eode. 


oi 


Finally, the variation of gamma that occurs across the 
contact surface in a wave rotor cannot be handled currently 


since E1DV2 assumes a single constant gamma throughout. 
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VI. CONCLUSIONS 


Towards the development of : one-dimensional code for 
wave rotor applications, methods for tracking and correcting 
conditions across a contact discontinuity, applying open and 
closed-end boundary conditions, accounting for shock wave 
and contact surface interaction were devised and were 
presented here in detail. The EULER1 Fortran Code [Ref. 2] 
was revised to become the E1DV2 code with the following 
additional capabilities: 

1) tracking of the contact surface and expansion wave 


2) imposing high pressure initially on the right side of 
the diaphragm 


3) correct jump conditions across the contact surface 
4) allow open boundary conditions with constant pressure 
specified and allow exiting of shock and expansion 


Waves 


5) allow closed boundary conditions and model shock and 
expansion wave reflections 


6) improved numerical accuracy from time zero to the 
maximum time step. 


The code was tested on the shock tube problem under four 
different initial and boundary conditions with excellent 
results. 
The following extensions are recommended to make the 
code suitable for wave rotor applications: 
1) Solve the Riemann problem at the moment when the shock 


and contact surface intersect 


SS) 


oe 


9) 


These 


Add additional code to track more than two discontinue 
ities and one expansion wave 


Incorporate a variable value of gamma into the code 


Update the characteristic curve approximation from 
linear to a higher order polynomial 


Add code necessary to handle open boundary conditions 
with inflow 


Improve on the numerical computation at time zero for 
low pressure ratios across the diaphragm 


Enable boundary conditions to be variable during 
program execution to allow wave rotor cycle design 


Compare computations using the code with experimental 


data where available 


Extend the code to two dimensiorce 


extensions may require various degrees or effort. 


However, the ability of the QAZ1D solution method to be 


extended rather simply to describe viscous multi-dimensional 


flow suggests that such efforts would be justified. 
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APPENDIX A 


DERIVATIONS OF EQUATIONS 


A. LIST OF VARIABLES 


Speed of sound 


Cp Specific heat at constant pressure 
Cyr Specific heat at constant volume 

= specific internal energy 

h Specific enthalpy 

P Static pressure 

Q Modified Riemann variable 


aR Reversible heat transferred 


Velocity magnitude 


R Modified Riemann variable 

Re Gas constant 

S Modified entropy (non-dimensional where S = S°Re¢) 
T Static temperature 

1S Time 

u Velocity relative to a standing shock wave 

Vv Specific volume 

W Incoming Mach Number relative to a stationary 


shock wave 


W Work 
p Density 
Y Ratio of specific heats 
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B. DERIVATION OF SHOCK JUMP EQUATION FOR HIGH PRESSURE ON 
THE RIGHT 


The analytical equation relating the non-dimensionalized 
extended Riemann variable, R, change through a normal shock 
is derived below similar to that for the Q variable change 
when high pressure is on the left [Ref. 2 :Appendix A]. 

Figure A.1 shows a shock moving to the left with 
velocity Vg. Subscript A is always associated with 
parameters to the left of any discontinuity, and B with 
those on the right. To enable normal shock relations to be 
used the system requires a velocity in the opposite direc- 
tion but of equal magnitude be imposed upon it. The 
coordinate system was defined such that vectors, such as 
velocity, directed to the right are vositive while those to 


the left are negative. 





Ursvedag 


Figure A.1 Shock Wave with High Pressure on Right 


GZ 


Since relative incoming Mach Number, w, relates the 
velocity downstream, or that region into which the sheck is 


traveling, to the speed of sound in that region [Ref. 10:pp. 


114-154] then 


(AT) 





w is a positive value because qa, though negative, is small- 
er in magnitude than the equally negative shock velocity, 
Vs. The speed of sound, A, is positive by definition. 

The appropriate extended Riemann variable in this case 


is R, defined as 


R = q- AS (A2) 


Since q is negative, this can be rewritten as 


ee acide 4S) ea) 


to emphasize that R is the Riemann variable associated with 
the lesser change across the shock [Ref. 3:pp. 18-21]. We 
adopt the pattern of non-dimensionalizing velocity by the 
speed of sound, pressure and density by corresponding values 
downstream of the shock. Using the entropy S, defined by 


Verhoff [Ref. l:p. 2] as 
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dp ~ 4a AB 
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AS A Ax B 
op ean @ it Pa Pay 
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2.2. b 7 eee 
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6 eI A I 

B “A B 2 B it lap ae 

= ——+ [1-,I I + ( ) a [In (s-(—) )] (A5) 
NA a ee Lees) Ss 

The ratios of pressure, density, and sonic velocity across a 


normal shock from Zucker [Ref. 7:p. 151] and Shapiro [Rea 


10:p. 118] in terms of Mach Number are 


SB ey 22s 2 ee (AG) 
Py v+1 y+l 
PR _ yt w* (A7) 
ms (y-Lw +2 
A 1/2 
i} oe i = iat Se ee A8 
a Tae ee -Dil  d ee (A8 ) 


The first term on the right side of equation (A5) can be ex- 


pressed in terms of w by applying simple continuity in mass. 
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Pa~uaAreay = 


Taking the areas as equal this becomes 


PpupAreap 
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By subtracting one from each side 
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Substitute equation (Al) for ug to get 


aoe | 2 (1 aoe) 


aay wang ee 
i cays 2(1 -w’) 
mn (y+l)w 
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(A9) 


(A10) 


(All) 


(A12) 


(A13) 


Combining equations (A6), (A7), (A8), and (A13) #£=into 


equation (45) gives 


Se AG 





4(Ul-wi) 2 2 lye 
22 7/2 
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C. DERIVATION OF CONTACT SURFACE JUMP EQUATION 
The analytical expression for the ratio of sonic 
velocity across a constant surface is derived, as follows: 


The first law of thermodynamics is stated as 


incremental work done ™ 
change of amount of on the 
internal energy = heat added - system by 

to system surrounding 


Ory, in diftierential form 


de = dqp t+ aw (A15) 


Internal energy is defined as 


then 
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dee ai —- bay — vap (A17) 


The incremental work is given by 


dW = pdv (A118) 


The heat addition is determined from the definition of modi- 


fied entropy, 


ds _1 SR 
pe) Y T 
dqr = - yrds (A19) 


Substituting equations (A17), (A18), and (A19) into equation 


(A15) yields 
dh - pdv - vdp = -yTdS - pav 
Canceling like terms, and rearranging gives 
dp = Yas + (2)an (A20) 
V V 
Enthalpy, h, is defined as 


Z (A21) 
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For a perfect gas the sonic velocity is given by 
A = (Rela 
thus 
oh Ee 5 Y Ret) 1/2 (art) (A22) 


Substituting equation (A21) and then (A22) into equation 


(A20), equation (A20) becomes 


y 


dp = irds + aiear, 
Wy 


aoe i 
seed 2 IG a Gord Keyes 


Magee a ab 2dA 
= vids + (=) Cpl ie (A23) 
For an ideal gas, 
t= Fence 
and 
Cy = Cy oe Rc 


Combining these two equations, and rearranging 
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enum Re (p (A24) 


Substituting equation (A224) into equation (A23), #£=and 
observing that pressure remains constant through the contact 
surface 

ZY 


ee.) eT GA 
— yds + F(R, Cay) a7 0 (A25) 


Eliminating T/v, equation (A25) becomes 


2 = ols (A26) 
Then S/R¢ becomes a non-dimensionalized entropy. Using the 
Metation of S = S/Re now for the non-dimensionalized forn, 
so that by integrating both sides 

B B 


ar J 8 = JF 


2 _ 
==z!S, -S,1 = ln A,/A3 


which can be written as 





Lr 
eA (A27) 
a pe Se 


Definitions used in this development, excluding modified 


entropy were from [Ref. ll:pp. 87-125]. 


209 


D. DERIVATION OF EXTENDED RIEMANN VARIABLE CHANGE ACROSS A 
CONTACT SURFACE 


An analytical expression for the change in the non- 
dimensionalized extended Riemann variable, Q, across a 
contact surface traveling right is derived here. Figure A.2 


depicts a contact surface moving right with velocity q. 


contact sumiace shock 
———— —— > 


IMoyngiajo| 





Figure A.2 Contact Surface Traveling Right 
Subscript Notation 


Values of parameters to the left of the interface are 
denoted by subscript A', and those to the right with sub- 
script B'. As illustrated, parameters to the left of the 
shock have subscript A, and to the right subscript B. All 
velocities are non-dimensionalized by sonic velocity 
immediately downstream of the discontinuity. Thus using the 


definition of the extended Riemann variable, 
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CPAs (A22) 


then 
Ons Q ' SINC +An Sa Gp! +Ap Sp 
Aas t Az 
Gar “An A ' 
A B A 
Se - A2g 
A a is Spi ( ) 
Since velocity is constant through a contact surface, SO 
that 


QO ' 0, ' A v 
A = 5s, -S,, 
B! Any A e 


Using equation (A27), this can be written 








Q, 170 0 _ ( 9 (Spi 7S, 1) 


a: 


] (A, /A,) (A31) 
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For a contact surface traveling to the left as 


illustrated in Fig. A.3, the derivation is entirely similar but 


shock contact surface 





: irl 
x 


Figure A.3 Contact Surface Traveling Left 
Subscript Notation 


using the R, extended Riemann variable with the result 


ees = Sexe (So Sony ies (A32) 
A Baa 
E. DERIVATION OF OPEN BOUNDARY CONDITIONS WITH CONSTANT 
PRESSURE 
The equation for density and temperature ratios in terms 


of pressure and entropy are derived below. The definition 


of modified entropy in non-dimensional form is 


2 1 P 
S = —- In (= A33 
5 ca 7 . 
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where P and » are non-dimensionalized ratics. Then 


2 


(S = Sy) (-y(y-1)) = In (P/p") 
ort ., 
= ay (ye 
y-1 mega 
Rearranging, = 
Ls --2p) (-y (y-))] 
oY = [Plex  ” 1] 
thus 


2 
(Cee ay ye)? <1/Y 
,  (S-Spyy-l) -y 
p = {Sep * )} (A34) 


using the perfect gas law:- 


(S 2p) (-¥ (y-1)) ae 


- 1 
oo (op lexp 
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Oe ie yee 
le) aa oT 


multiply each side by T and divide by 97” then 
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APPENDIX B 


OPERATION OF E1DV2 ON THE NPS VM/CMS_ SYSTEM 


A. TERMINAL REQUIREMENTS 

Terminal requirements depend on the desired output. Any 
terminal, such as the IBM 3278, connected into the VM/CMS 
system is adequate for a tabular listing or Disspla Metafile 
Se data. The Disspla Metafile stores graphical data for 
display using DISSPOP commands. If output of the plots ona 
monitor screen is desired, an IBM 3277-TEK618 dual screen 
terminal must be used. Table B.I lists terminal require- 


ments for different outputs. 


TABLE B.I 
TERMINAL AND OUTPUT 
OUTPUT TERMINAL 
Tabular listing of data Any terminal 


Graphical data in format 
Bene VERSATEC printing Any terminal 


Graphical data plotted 
on screen ieMeSe27 7-TEKG1S 


B. PROCEDURE 


1. Editing Variables and Program Setup 


To change the initial and boundary conditions, 


graphical output, and grid size to run the program under 
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different conditions, the following must be done. Five 
pick the appropriate terminal from Table B.I and log on. 


Enter the editor by issuing the command 
xX E1DV2 FORTRAN A 


Table B.II lists the code variables that will require 
editing and their meaning. In addition, to these variables 
if graphical plots are to be outputed then enter the 
"BORDER" and "EXACT" subroutines and edit the lines: 
MrERST ORDER N = 222" 
"DENSITY RATIO = 7??? TEMP RATIO = ???" 
"PRESSURE RATIO = ??°?" 
To issue the output graphs to the screen comment out the 
lines "COMPRS", and use 
CALL TEK618 
Otherwise, if a Disspla Metafile of the graphics plot is 
desired, comment out the "TEK618" line, and use 
CALL COMPRS 
These lines are in the "Set up graphics plot of variable" 
loop in the main program. 

To store these changes, hit the enter key, and type 
"file". Select the enter key once more. The program is now 
ready to be compiled and executed. 

2. Commands 

To compile and execute the E1DV2 code on the VM/CMS 


system perform the following commands exactly as typed: 


J 


TABLE B.IfT 


EDITABLE VARIABLES 


Variable Bdit 
DIMENSION statement set the dimension of the arrays 
for arrays A(N),.-.., equal to the number of nodes 
XARRAY (N) in the grid 
N set to the number of nodes in the 

grid 
GRAPHS set to: 0O for tabular listing of 
data 


ites plot of density, 
entropy, pressure and 
velocity distributions 

2 for plot of exact density 
distribution compared to 
computed density 


distribution 
SKIP set to number of time steps between 
Calls to output routines 
JSTOP set to maximum number of time steps 
TRE set to initial temperature ratio 


across diaphragm 


PRI set to initial pressure ratio across 
diaphragm 

DRI set to initial density ratio across 
diaphragm 

LL set to initial velocity left of 
diaphragm 

QRI Scoeerbouinihctal velocity right of 
diaphragm 

LBDPRI set to initial pressure ratio 


across left boundary 


LBDTRI set to initial temperature ratio 
across left boundary 


LBDORI set to initial density ratio across 
ijeft boundary 


ay 


RBDDRI 


RBDPRI 


RSDIRE 


EE 


LWPRES 


LBNDRY 
RBNDRY 


LBDPRS 
RBDPRS 


XINIT 


VHEAD 


VTAIL 


VCDE 


VSE 


DLCD 


DLSH 


TABLE B.IL (CONTINUED) 


set to initial density ratio across 
right boundary 


set to initial pressure ratio across 
right boundary 


set to initial temperature ratio 
across right boundary 


set to desired value of gamma 


set to desired error tolerance for 
calculation of characteristic slope 
(l1.e., 0.1D-8) 


set to: 2 for low pressure on right 
side of diaphragm 
3 for low pressure on left 
side 


set to: 1 for closed boundary 
O for open boundary 


set to: 0 for constant pressure at 
left boundary 
1 for pressure that adjusts 
at the left boundary 
set to location of diaphragm 


set to exact velocity for head of 
expansion wave 


set to exact velocity for tail of 
expansion wave 


set to exact velocity for contact 
surface 


set to exact velocity for shock 
wave 


set to exact density behind contact 
surface 


set to exact density behind shock 


ES 


TABLE B.II (CONTINUED) 


SIGMA (1,2) set to diaphragm location (i.e., 
SIGMA (2,2) 0.5D00) 

SIGMA (3,2) 

SIGMA (4,2) 

Y There are four cases where Y 


appears in the program 

edit as follows: 

Comment out Y = (Integer#), use: 
first, Y = (N+1)/2 if LWPRES = 2 
second, Y = (N+3)/2 

and 

third, Y = (N-1)/2 if LWPRES = 3 
fourEneee ean t+ 1) /2 

for diaphragm at 0.5 

Comment out those above, if dia- 
phragm at another node. Set 


first, Y = (Integer #) of node for 
LWPRES = 2 

second, Y = (# + 1) 

and 


third, Y = (Integer # - 1) of node 
for LWPRES = 3 
fourth, Y = (#) 


JEL, 


1) Increase the virtual memory by entering 
DEFINE STORAGE 1M 
2) To return to CMS environment enter 
I CMS 
3) To compile the program enter 
FORTVS E1DV2 


The screen will display messages as it compiles each 
routine and when finished a ready symbol appears. 


4) To execute the program, enter 
DISSPLA E1DV2 
The message 


--- YOUR FORTRAN PROGRAM IS NOW BEING LOADED... 
.- EXECUTION WILL SOON FOLLOW... 


should appear, followed by 
we ADCULION BEGINS 

If at a TEK618 terminal with GRAPHS equal to 1 or 2 
then the screen on the TEK618 will begin plotting the 
selected graph. A 

-»-press ENTER €0O continue 
message will appear on the 3277 terminal. If a copy of the 
plot is desired, do so now before pressing the enter key. 
After pressing the enter key on the 3277 terminal, the plot 
will be erased and the program will terminate. Proper 
termination will result in 
END OF DISSPLAY 9.2 #77 VECTORS IN] treet 


appearing, followed by a ready message. 
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If GRAPHS was set to 0, then proper terminaticn 
would be a ready message. The tabular listing of pressure, 
density, velocity, entropy, and Riemann variables will be in 
"FTLE FTO9FOO1." The exact location of the shock, contact 
surface, and expansion wave with elapsed time will be in 
wetie FTOSFOOl." The computed location of the shock, 
contact surface, and expansion wave with elapsed time will 
Benin “FILE FTI10FO0O1." 

A second method to compile and execute the program, 
plus provide the files with a name is to create the follow- 
ing EXEC file on the user's disk. 

foe 2 DISK FILEOO LISTING A( PERM 
Poe LO DISK FILELO LISTING A( PERM 
FI 8 DISK FILEO8 LISTING A( PERM 
FORTVS &1 
A possible file name and required file type would be 
RUN EXEC 
To compile the program and define the output files, enter 
RUN ELDV2 
After compiling is finished, and the ready message appears, 
enter 


DISSPLA E1DV2 


The program executes as outlined before. 
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APPENDIX C 


FLOWCHARTS 


Flowcharts for major routines in E1DV2 are shown. 
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Figure C.1 Main Program Flowchart 
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Figure C.7 (Contamued) 
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Figure C.2 "DBURST" Subroutine Flowchart 
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Figure C.2 (Continued) 
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Figure C.3 "TRAK" Subroutine Flowchart . 
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Figure C.3 (Continued) 
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Figure C.3 (Continued) 
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Figure C.3 (Continued) 
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Figure C.4 "SWEEP" Subroutine Flowchart 
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Figure C.4 (Continued) 
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Figure C.4. (Continued) 
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Figure C.5 General "COND1, 2, 3, and 5" Subroutine Flowchart 


Input Data 











no yes 


ode 
oni Be Mic 





set : set 
Rnode(1) =I Lnode(1) =! 
Rnode(2) = snock | Lnode(2) = snock 


Rnode(S) = cntact | Lnode(3S) = cntact 
Rnode(4)=j | Lnode(4) = | 





Figure C.6 "“COND4" Subroutine Flowchart 
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Figure C.7 "CORRCT" Subroutine Flowchart 
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Figure C.7 (Continued) 
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Figure C.7 (Continued) 
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Figure C.8 "BONDRY" Subroutine Flowchart 
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Figure C.8 (Continued) 
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Figure C.9 "SRFLCT" Subroutine Flowchart 
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ANNANDANAAAAAAANAAANAAAAAAAANAADAADAAAAAAAAADAAANAANAAAAADAAAAAAADHHAAAA 


APPENDIX D 


E1DV2 FORTRAN LISTING 


Tob tototaketotaketetetetotetetetataheiahatatatetetehahetetetetatetetetetetateteteteateteletete! 


EULER-10 
VERSION 2 
{E1D0V2) 


THIS PROGRAM SOLVES THE EULER EQUATIONS 
EXPRESSED IN A QUAZT-ONE DIMENSIONAL 
STREAMLINE COORDINATE SYSTEM. 


AUTHOR - LT. D.T. JOHNSTON,FEB 1987 


BASED ON THE EULER1 CODE BY 
T.F. SALACKA, DEC 1985 


FEATURES OF THIS VERSION (2) 

¢ ORDER OF SPATIAL DERIVITIVES © FIRST 
# NUMBER OF SPATIAL DIMENSIONS - ONE 

# DISCONTINUITIES TREATED: 


a a a a a a a 2 o> o> oe ae ae ae ae oe Se Se Se oe Se oe ae ae oe 
a a oe oe oe oe oe oe oe oe oe Se Se ae oe Se oe oe Se oe oe 


SHOCKS e YES 

CONTACT OISCONTINUITIES - YES 

EXPANSION WAVES e YES 
@ HIGH PRESSURE SIDE 

LEFT e YES 

RIGHT © YES 


DE OG 3 SEE FG TG OE FG FG 2G at 2G FE 2G 0 96 OG OE OE OG FE OG 90 9G 20 FG OG FE Ot BE DE-DE 00 30 BE 26 26 OE 26 26 06 20 3G 2 38 30 Ot 


SHHHHEHEEHEFEESEHEHSEEHHESEEHEDEFESHES SHEESH E HEHE SHEE DHS 


+ 
+ CONVENTIONS AND OEFINITIONS + 
« ¢ 


HHH G HFG 4H HHH HHH S444 44444444449 HHHHH HES 
eoe--= NON-DIMENSIONING CONVENTION----------- 


ALL VELOCITIES NON-OIM. BY THE SOUND SPEED ON 
THE LOW PRESSURE SIOE OF THE OTAPHRAGM. 

ALL PRESSURES, DENSITIES, AND TEMPERATURES ARE 
NON-OTHM. 8Y THEIR INITIAL VALUES ON THE LOW 
PRESSURE SIDE OF THE OIAPHRAGM, 

SPACITAL DISTANCE IS NON-DIM.BY OVERALL LENGTH. 

ENTROPY IS NON-OIM. BY THE GAS CONSTANT, R. 

TIME IS NON-DIM. BY {LENGTH/SOUND SPEED). 

VELOCITIES AND DISTANCES ARE DEFINED POSITIVE 
TO THE RIGHT. 


eerese-s= © SUBSCRIPT NOTATION------------- o--- 


- SPACIAL NODE (1 TON FROM LEFT TO RIGHT) 
TIME LEVEL (0 IS THE INITIAL CONDITION) 
DENOTES WHICH CHARACTERISTIC WAVE IS BEING 
OEALT WITH: 

1 @ Q+A 22 38 Q-A 
L = DENOTES WHICH TYPE OF DISCONTINUITY IS 
BEING DEALT WITH: 

1 = SHOCK 2 ® CONTACT DISCONTINUITY 

3 2 HEAD OF EXPANSION HAVE 

4 = TAIL OF EXPANSION HAVE 


A Go 
é 
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A - SPEED OF SOUND 
%A # DENOTES THE VALUE OF A VARIABLE AT THE NODE 

TO THE LEFT OF A DISCONTINUITY. % CAN 

BE ANY VARIABLE NAME. 
AR = THE RATIO OF SOUND SPEED ACROSS A 

SHOCK, A/BI SHOCK MOVING RIGHT ),B/A(SHOCK MOVING LEFT) 
*B - DENOTES THE VALUE OF A VARIABLE AT THE NODE 

TO THE RIGHT OF A DISCONTINUITY. 
BDRY - 3 DENOTES LEFT BOUNDARY,2 THE RIGHT BOUNDARY 
COUNT <= COUNTER FOR GRAPHICS ROUTINES 
DARRAY =- ARRAY OF DENSITY FOR PLOTTING 
DELT - TIME STEP 
DENSITY BEHIND THE CONTACT 
DISCONTINUITY IN THE EXACT SOLUTION. 
DENSITY BEHIND THE SHOCK IN THE 
EXACT SOLUTION. 
DQ - THE JUMP IN VELOCITY ACROSS THE SHOCK 

DIVIDED BY THE SOUND SPEED AT B(RIGHT) OR A(LEFT) 
ORI - INITIAL DENSITY RATIO ACROSS THE SHOCK 
DESIRED PRECISION FOR CHARACTERISTIC CALCULATIONS 
G - GAMMA (RATIO OF SPECIFIC HEATS) 
GRAPHS = FOR GRAPHICAL OUTPUT, O=NONE (TABULAR) 

1=PLOTS ALL VARIABLES 
2=COMPARES DENSITY WITH EXACT SOLUTION 


oO 
Ca 
(2) 
Oo 
( 


Oo 
tr 
W 
=x 
( 


m 
m 
Q 


Gl - 1/(G-1) 
G2 - 2/(G-1) 
H - 1/(N-1) 
HALT = TERMINATES PROGRAM IF 1,SET BY CONDITIONS NOT CODED 
id - NUMBER OF THE NODE TO THE RIGHT OF A 
DISCONTINUITY. 


JSTOP - NUMBER OF TIME LEVELS TO BE CALCULATED 

LBDDR - LEFT BOUNDARY DENSITY RATIO 

LBODRI = LEFT BOUNDARY DENSITY RATIO AT TIME ZERO 

LBDOPR = LEFT BOUNDARY PRESSURE RATIO 

LBDPRI - LEFT BOUNDARY PRESSURE RATIO AT TIME ZERO 

LBDPRS - VALUE OF O DENOTES CONSTANT PRESSURE AT LEFT BOUNDARY 
»1 DENOTES ADJUSTABLE PRESSURE AT THE LEFT BOUNDARY 

LBOTR = LEFT BOUNDARY TEMPERATURE RATIO 

LBOTRI - LEFT BOUNDARY TEMPERATURE RATIO AT TIME ZERO 

LBNORY - DENOTES LEFT BOUNDARY CONDITION, OPEN OR CLOSED 

LNODE - ARRAY OF LEFT MOST NODE TO BE CORRECTED IN CORRCT 

LWPRES =- DENOTES WHICH SIDE OF DIAPHRAGM HAS LOW PRESSURE 

N - NUMBER OF SPACIAL NODES (ODD NUMBER) 

ND - DOUBLE PRECISION VALUE OF N 

NEW*(T J- STORED VALUES OF ** FOR THE NEXT TIME LEVEL 

PARRAY - ARRAY OF PRESSURES FOR PLOTTING 

PRI - INITIAL PRESSURE RATIO ACROSS THE SHOCK 

PLTCNT COUNTER FOR GRAPHICS ROUTINES 

Q ABSOLUTE FLUID VELOCITY 

QARRAY - ARRAY OF VELOCITIES FOR PLOTTING 

QLBD - INITIAL VELOCITY AT LEFT BOUNDARY 


QLI - INITIAL VELOCITY LEFT OF THE DIAPHRAGM 
QRBD - INITIAL VELOCITY AT RIGHT BOUNDARY 
QRI - INITIAL VELOCITY RIGHT OF THE DIAPHRAGM 
QQ - Q+A¥S (EXTENDED RIEMANN VARIABLE ) 


RBDDR - RIGHT BOUNDARY DENSITY RATIO 

RBODORI - INITIAL RIGHT BOUNDARY DENSITY RATIO 

RBDPR - RIGHT BOUNDARY PRESSURE RATIO 

RBDPRI - INITIAL RIGHT BOUNDARY PRESSURE RATIO 

RBDPRS - VALUE OF O DENOTES A CONSTANT PRESSURE AT RIGHT 
BOUNDARY, WHILE 1 IS FOR ADJUSTABLE PRESSURE 

RBDTR =- RIGHT BOUNDARY TEMPERATURE RATIO 

RBOTRI - INITIAL RIGHT BOUNDARY TEMPERATURE RATIO 

RBNDRY = DENOTES RIGHT BOUNDARY OPEN OR CLOSED 

RNODE - ARRAY FOR RIGHT MOST NODE TO BE CORRECTED IN CORRCT 

RR - Q-A®S CEXTENDED RIEMANN VARIABLE ) 

S - ENTROPY : 

SARRAY - ARRAY OF ENTROPY FOR PLOTTING 
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SIGMA - SPATIAL LOCATICN OF DXSCCNTINUITIES 
SIGMA(L,J) WHERE L INCICATES THE TYPE OF 
DISCONTINUITY AND J INDICATES THE 
TIME LEVEL$ 1 - CURRENT LEVEL 
2 - LEVEL BEING CALCULATED 
SK =- INTEGER THAT DENOTES RELATIVE LOCATION OF SHOCK NEAR 


BOUNDARIES 
SKIP - VARIABLE WHICH INDICATES HOW MANY TIME STEPS BETWEEN 
CALLS TO OUTPUT ROUTINES 
T - TIME SINCE INITIAL CONDITIONS 
TRI - INITIAL TEMP. RATIO ACROSS THE SHOCK 


VHEAD = VELOCITY OF THE HEAD OF THE EXPANSION 
WAVE FOR THE EXACT SOLUTION. 

VTAIL - VELOCITY OF THE TAIL OF THE EXPANSION 
WAVE FOR THE EXACT SOLUTION. 


VCDE - VELOCTIY OF THE CONTACT DISCONTINUITY 
FOR THE EXACT SOLUTION. 
VS - THE SHOCK SPEED( POSITIVE RIGHT, NEGATIVE LEFT) 
VSE - VELOCITY OF THE SHOCK FOR THE EXACT 
SOLUTION. 
W - MACH NO. RELATIVE TO A STANDING SHOCK 


XARRAY = ARRAY OF SPATIAL POSITIONS FOR PLOTTING 
XEXACT = ARRAY OF SIX X VALUES FOR THE EXACT SOLUTION. 
XINIT = INITIAL POSITION OF DISCONTINUITY FOR 

EXACT SOLUTION PLOTTING. 


X2 - LOCATION OF NODE TO RIGHT OF DISCONT. 
ALONG THE SPACTAL AXIS. 
Y - (N+1)/2 


YEXACT - ARRAY OF SIX DENSITY VALUES FOR THE EXACT SOLUTION. 


30% OTHER VARIABLES ARE DEFINED IN THE 
SUBROUTINES WHERE THEY ARE USED 33% 


SHESHHSESESHEHSES HEPES SESE SESESE EHS ESE SESESEH EHH HEE HEE EHH 


+ ~ 
~ PROBLEM SET - UP - 
+ + 


FEEEHEEEEEEEEESE SEE ESEE SESE FES EE EEEFESES HEHEHE HHH 
THE PARTICULAR PROBLEM FOR THIS VERSION IS: 


SHOCK TUBE, SINGLE CENTERED DIAPHRAGM WITH 
HIGH PRESSURE SIDE TO THE RIGHT. 


BOUNDARY CONDITIONS = LEFT END CLOSED,RIGHT END OPEN 
--------- VARIABLE DECLARATIONS <------------- 


DIMENSION 12(4),X2(4),XEXACT( 6) »YEXACT( 6 ),LNODE(4),RNODE(4), 
C SIGMA(4,2) 


+++++++ USER INPUT REQUIRED HERE +44+44+44444444 


DIMENSION A(101),Q(101),QQ(101),RR(101),S(101); 
Cc NEWRR( 101 ),NEWS(101),NENQQ( 101), 
C PARRAY(101),DARRAY(101),SARRAY(101), 
C QARRAY(101),XARRAY(101 ) 


INTEGER I>J»N»JSTOP >Y,GRAPHS »COUNT ,PLTCNT »BORY »SK,RBNDRY » 
SKIP ,IT2,LWPRES >HALT ,LNODE ,RNODE ,LBNDRY ,LBDPRS»RBDPRS 

DOUBLE PRECISION TRI>,PRI,QLI,QRI,ORI,G,G1,G62,SIGMA,EE ,NENQQ, 
DELT »H»ND,X2,AR,W,D0Q,VS>T »A»Q,QQ>,RR»S,»NEWRR > NEWS» 
LBOPRI,LBOTRI,LBDORI,LBOPR,LBDDR,LBOTR,QLBD, 
RRA, QQA.,AA,SA,QA,RRB ,QQB,AB>SB>GB, 
RBOPRI,RBOTRI,RBOORI,RBOPR,RBOODR ,RBOTR,QRBD, 
QCS ,VTEW,VHEW 

REAL VTAIL,VCDE,VSE ,DLSH,DLCD ,»XINIT ,VHEAD, 


oO 


aaaan 


145 


Cc XE XACT » YEXACT ,PARRAY » DARRAY »SARRAY »QARRAY »XARRAY 
COMMCN AR,DQ,VS>W 


Cc 
C £+---- ENTER THE APPROPRIATE VALUES BELOW ------ 
Cc 
N= 101 
GRAPHS=1 
c , 
C @---- FOR GRAPH = 1 OR 2 MUST ENTER CHANGES IN SUBROUTINE ----- 
C #£ e222 BORDER AND SUBROUTINE EXACT ----- 
C 8 £+e2-- LINES “FIRST ORDER No = 777" 
C £#-e2- “DENSITY RATIO = 777 TEMP RATIO = 777" 
C ofe-+2 “PRESSURE RATIO = ?77?°* 
C 
SKI P=18 
JSTOP=101 
TRI=1.0D00 
PRI=5.0000 
ORI=5.0D00 
QLI=0.0000 
QRI=0.0D00 
LBDPRI = 1.000 
LBDTRI = 1.000 
LBODRI = 1.D00 
RBODRI = 1.000 
RBDPRI = 4.BD00/5.D00 
RBDTRI = 1.000 
G=1.4D00 
EE=0.1D-8 
c 
C e---- DENOTE LOW PRESSURE SIDE BY SETTING ----- 
C ----- LWPRES = 2 IF LOW PRESSURE ON RIGHT 
C errr" LWPRES = 3 IF LOW PRESSURE ON LEFT 
c 
LWPRES = 3 
c 
C -c--- SET BOUNDARY CONDITIONS BY SPECIFING OPEN OR CLOSED ----- 
C ----- LBNORY: OPEN = 0, CLOSED = 1 (FOR LEFT BOUNDARY) ----- 
C ----- IF OPEN SPECIFY IF PRESSURE IS TO BE MAINTAINED AT LBDPRI ----- 
C --- OR IF IT CAN ADJUST TO PREVENT ANY WAVES FORMING AT THE BOUNDARY - 
C rere LBDPRS: CONSTANT = 0, ADJUSTABLE = 1 ----- 
C e---- DO THE SAME FOR RIGHT BOUNDARY;3 RBNORY,RBOPRS ----- 
Cc 
LBNDRY = 1 
LBOPRS = 1 
RBNDRY = 0 
RBDPRS = D0 
C = _ eeeerrrre EXACT SOLUTION VALUES ---e-------- 
XINIT=0.50 
VHEAD= 1.0 
VTAIL= 0.310557 
VCDE=-.574487 
VSE=-1.402346 
OLCD=2.713115 
DLSH=1.69344 
SIGMA(1,2)}=0.50000001000 
SIGMA(2,2)=0.50000001000 
SIGMA(3,2)=0.50000001000 
SIGHA(4¢,2)=0.50000001000 
Cc 
Cc 4444444444 END OF USER INPUT AREA 44444444 
C 
SK = 0 
T=0.0D00 
DO 10 I=1,4 
I2( TI }=0 
LNODE(I} = 0 
RNODE(I) = O 


X2(1I)=0.0060 
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SIGMA(T,1)=D0.0DC0 
10 CONTINUE 
ND=DBLE(N) 
H=1.D0D0/(ND-1.060) 
DO 11 I=1,N 
A(I )=0.0D00 
Qi I )=0.0000 
NEWHQQ(I )=0.0000 
NEWRR( I )=0.0000 
NEWS( TI )=0.0000 
QARRAY(I )=0.0 
PARRAY(I )=0.0 
DARRAY(I )=0.0 
SARRAY( I )=0.0 
XARRAY( TI }=FLOAT( I-1 )¥SNGL(H) 
11 CONTINUE 
DELT=2.0D00 
Gc 
C w----LOAD INITIAL REIMAN VALUES INTO NODE LOCATIONS,FIRST FROM NODE--- 
C ----- 1 TO MIDPOINT (Y), AND THEN FROM Y TO N. NOTE IF SHOCK DOES NOT-- 


AR=1.0000 
DQ=0.0000 

W=1.DD00 

VS=0.0000 
QCS=0.D00 
VHEW=0.D00 
VTEW=0.000 
G1=1.000/(G-1.000) 
G2=2.D000/(G-1.D00 ) 


C ----- FLOWN RIGHT ----- 


IF(LWPRES.EQ.2) THEN 
LBDDOR DRI * LBDDRI 
LBDPR PRI * LBOPRI 
LBOTR TRI * LBOTRI 
QLBD = QLI 

= RBODRI 
RBOPR = RBDPRI 

= RBDTRI 
QRBD = QRI 

Y = (N+1)72 


DO 12 I=1,Y 
S(T }=G2-(G61/G )*DLOG( PRI/( (ORI )**G) ) 
QQ( TI J=QLI+DSQRT( TRI I¥*S(T) 
RR(I }=QLI-OSQRT( TRI )¥S(T) 
12 CONT INUE 
Y=(N4#3)/2 
Cc Y=39 
DO 13 I=Y,5N 
S$(I)=G2 
QQ( I )=QRI+S(T ) 
RR(IV=QRI-S(T) 
13 CONTINUE 
ELSE 
C ----- FLOW LEFT ----- 
LBDOOR LBDDRI 
LBDOPR LBOPRI 
LBDTR LBDTRI 
QLBD = QLI 
RBDDR RBDDRI * DRI 
RBDPR RBDPRI * PRI 
RBOTR RBOTRI * TRI 


147 


QRBD = QRI 


Y = (N-1)72 

Y=13 
00 14 I=1,Y 
$(I)=G2 


QQ(T J=QLI+S(T) 
RR(I)=QLI-S(T) 
14 CONTINUE 


Y = (N#1)/2 
Y= 14 
90 15 I=Y,N 


S(I )=G2-(G1/G )*¥OLOG(PRI/( (ORI )**G )) 
QQ( TIT J=QRI+OSQRT( TRI )*S(T) 
RR( I J=QRI-OSQRT( TRI I*S(T ) 
15 CONTINUE 
END IF 


IF (GRAPHS.GT.0O) THEN 
CALL COMPRS 

CALL TEK618 
CALL HWROT( ‘AUTO‘ ) 
CALL HWSCAL( ‘SCREEN’ ) 
IF (GRAPHS.EQ.2) THEN 

CALL PAGE(11.0,8.5) 
ELSE 
CALL PAGE(8.5,11.0) 

ENO IF 
IF (GRAPHS.EQ.1) THEN 
CALL BORDER( JSTOP ) 
ENO IF 

ENO IF 

HALT = 0 

J=1 

COUNT =1 

IF (GRAPHS.EQ.1)} THEN 
CALL PLOT(J,JSTOP »N,QQ,RR»S»H»XARRAY »PARRAY » 

#OARRAY »QARRAY »SARRAY »G,G1,G2) 
ENO IF 


~---- BURST DIAPHRAGM ----- 


CALL TIME(N,QQ,RR,S,OELT,H) 
CALL OBURST(N,H,QQ,RR,S,G,G1,G2,D0ELT,I2,X2,W,AR,0Q,VS,LWPRES, 
Cc SIGMA,A,Q} 


IF (GRAPHS.EQ.0) THEN 
CALL LIST(N,SIGMA,QQ,RR>S,G,G1,62,J,T DELT,VS»QCS,VTEW,VHEN, 
C RINIT »VHEAD ,»VTAIL,VCOE ,VSE ) 
ENO IF 


<== BEGIN CALCULATION FOR JUMP TO NEXT TIME ANO CONTINUE -<---- 
=--=— UNTIL EITHER JSTOP REACHED OR SHOCK MEETS CONTACT SURFACE ~---- 


16 =IF (J.EQ.JSTOP) GOTO 18 
PLTCNT=J/SKIP 


CALL TIME(N,QQ,RR»S,OELT,H) 


CALL TRAK(N,SIGMA,H,QQ,RR»S,G,G1,G62,0ELT,I2,X2»N,AR,OQ,VS>J> 
Cc LWPRES,QCS,VHEW, VTEW ) 


CALL SWEEP(N,H»SIGMA,QQ,RR,»,S,OELT,EE »Q>A ,NEWQQ ,NEWRR »NEWS »I2,G62, 
J,LNOOE ,RNOOE ,HALT,LBNORY,LBOPRS,LBOPR,LBOTR,LBOOR, 
QLBO ,G,G61,RBNORY »RBOPRS,RBOPR,RBOTR,RBOOR ,QRBO >SODRY » 
SK ) 


aan 


IF(LNOOE(4%).LT.J) THEN 
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LNOCE(1) «= 
RNODE(1) = 
GO TO 17 


0 

0 

END IF 
IF(RNODE(G).LT.LNODE(4)) THEN 


RNODE(1) = 0 
END IF 


CALL CORRCT( LNODE ,RNODE »5N,SIGMA,H,QQ,RR,S,G,G1,62,I2,X2,N,AR,DQ, 
Cc VS »A,Q) 


17 IF (SIGMA(1,2).GE.1.000) THEN 
IFC RBNORY.EQ.0) THEN 
IF(SIGMA(1,2).NE.3.000) THEN 
SK = 2 
CALL BONDRY(Q(N),Q(N-1),QRBD,A(N),A(N-1),QQ(N),QQ(N-1), 
RR(N),RR(N-1),S0N),S(N-1),H,EE,DELT, 
RBNORY ,RBDPRS ,»RBOPR,RBOOR ,RBDTR>J» NEWQQ(IN), 
NEWRR(N),NEWS(N),G,G1,G2,HALT,B8DRY ,SK ) 


eHaAn 


END IF 
ELSE 
CALL EXTRAP(RRIN=-1),RRIN-2),QQUN=-1),QQUN-2),S(N-1), 
Cc S(N=2),H»H»RRA»QQA SA >,AA QA } 
CALL SRFLCT(QQA,RRA»,SA,SIGMA,VS ,DELT,LWPRES, 
Cc RR(N),QQ(N),S(N),QUN),A(N},G6,G61,G62) 
END IF 
END IF 


IF (SIGMA(1,2).LE.0.900) THEN 
IF(LBNORY .EQ.0} THEN 

IF(STIGMA(1,2).NE.-2.000) THEN 

BORY = 3 

SK=2 

CALL BONDRY(Q(1),Q'2),QLBD,A(1),A(2),QQ(1),QQi2), 
RR(1),RRO2),501),502),H,EE,OELT > 
LBNORY,LBOPRS,LBOPR,LBODR,LBDOTR>»J,NENQQI 1), 
NEWRR(1),NEWS(1),>G,G1,G2,HALT ,BDRY,SK ) 


aan 


END IF 
ELSE 
CALL EXTRAP(RRE2),RRE3),QQI2),QQ(3),S(2), 
Cc S(€3),5H»,H,RRB,QQB,SB,AB,GB } 
CALL SRFLCT(QQB,RRB>SB,SIGMA,VS,DELT,LWPRES, 
Cc RRO1),QQI1),501),Q01),A01),6,61,62) 
END IF 
END IF 
T=T+DELT 
Cc 
C =---- OUTPUT DATA ----- 
c 
IF ((COUNT.EQ.PLTCNT*SKIP ).ANO.(GRAPHS.EQ.1)) 
Cc THEN 
Cc IF((J.GT.55)) THEN 
CALL PLOT( J» JSTOP ,N»,QQ>RR»S>H»XARRAY »PARRAY » 
CDARRAY ,QARRAY ,SARRAY »G»G1 5,62) 
Cc END IF 
END IF 
IF ( (COUNT .EQ.PLTCNT*SKIP).AND. (GRAPHS .EQ.0)) 
Cc THEN 
CALL LIST(N,SIGMA,QQ,RR»S,G,G1,62,J,T ,DELT»VS,QCS,VTEW,VHEW, 
Cc XINIT,VHEAD,VTAIL,VCOE »VSE ) 
END IF 
IF ( (COUNT .EQ.PLTCNT*SKIP).ANOD.(GRAPHS.EQ.2)) 
Cc THEN 
IF (J.GT.50) THEN 
CALL EXACT(N,XINIT,T,VHEAD ,VTAIL »VCDE »VSE ,DLCO ,DLSH,QQ,RR>S>H, 
CXARRAY ,DARRAY »G>»G1,G2>,0RT ) 
ENO IF 
END IF 
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TFCHALT.EQ.1) GO TO 18 

J=J+1 

COUNT =CCUNT +1 

GO TO 16 

18 CALL DONEPL 

END 
3300000000800 BBB BB BRB ERERE 
336-36 DEDEDE DE DEE DEE 36 3 IE 3 IE DIE 9 3 36 EE EE DEE IE 3 9 EE EE 9 DEE 96 IEE IEEE 23 
HEE 3 IE DE IE DEE 3 IE IE SESE 3 ESE IE 3 EE SEE DEE 96 EEE SE IESE DE 96 36 36 3 IE DOSE EE 3 3E 96 98 3 


SUBROUTINE LIST(N,SIGMA ,QQ,RR»S,G,G1,G62>,J,T »DELT»VS,QCS,VTEW, 


C VHEW, XINIT »VHEAD »VTAIL ,VCDE ,VSE ) 
FHEFHEHEEE HEHE FEE EE ESES ES ES EFESE ESE SESES HE HHH HHH tet 
* $ 
+ TABULAR RESULTS SUBROUTINE + 
+ + 


PHESH HHH HHS $$$ 4444444464444 4444444446444444444444 
w---------- VARIABLE DEFINITIONS --------------- 


DENS - DENSITY 
PRESS = PRESSURE 
TEMP =- TEMPERATURE 


INTEGER I>,J>N>L 

DIMENSION SIGMA(4,2 ),QQ(N),RRIN),SON) 

DOUBLE PRECISION SIGMA,QQ,RR»S,PRESS»VTEW,VHEW>QCS » TENXE » 
Cc TEMP ,DENS,G,G1,G2,Q,T ,OELT,VS,HEWXE , 
C XINIT ,VHEAD > VTAIL, VCDE » VSE »SKXE »CSXE 


WRITE(9,%) "TIME LEVEL’>J>’ ELAPSED TIME IS',>T 
WRITE(9,%) ‘TIME STEP IS',DELT,’ SHOCK VELOCITY IS',VS 
WRITE(9,%#) ‘CONTACT SURFACE VELOCITY IS',QCS 
WRITE(9,4) ‘HEAD EXPANSION WAVE VELOCITY IS',VHEW 
WRITE(9,#) ‘TAIL EXPANSION WAVE VELOCITY IS°*,VTEW 
WRITE(9,#) ° ' 
WRITE(9,%) ° NODE VELOCITY DENSITY 
CPRESSURE ' 

WRITE(9,#) ° ° 
DO 61 I=1,N 
TEMP=(QQ( TI )-RR(I ))*(QQ(T)-RROI)D)/7(4.000*S(I)*S(T)) 
DENS=((1.000/TEMP )¥DEXP(G*(1.D00-G )*(S(I )-G2)) )**(-G1) 
PRESS=TEMP*QENS 
Q=(QQ( I )+RRII))/2.0000 
WRITE (9,65) I,Q,O0ENS,»PRESS 

65 FORMAT (4¢X,I2)7X,F12.6,7X,F12.6)7X»yF12.6) 

61 CONTINUE 
WRITE(9,#} ° ° 
WRITE( 9,4) ' ' 
WRITE(9,%) ' NODE QQ RR 

CMODIFIED S° 

WRITE(9,#) ° ° 
DO 62 I=1,N 
WRITE (9,66) I,QQ(I),RRII),S(T) 

66 FORMAT (4X%,12,7X,F12.6»)7X»F12.6)7X,F12.6) 

62 CONTINUE 
WRITE(9,%) ' 
WRITE(9,%#) ° 
WRITE(9,*) ° DISCONTINUITY LOCATIONS AT TIME LEVEL*,»J 


WRITE(9,*) 
WRITE(9,#) 
DO 63 L=1,4 
WRITE(9,%) L,’ * SSIGMA(L,2) 
63 CONTINUE 
WRITE(9,*)} * ' 
WRITE(9,%) '-------------------- wo 2---2-------------------------- 


Us = LOCATION’ 


WRITE(9,*) ° ° 


qaeaoanaanaanaannanan © 


aAeagandgNnanannanannanann 


WRITE(9,%) ° ° 
IF(J.EQ.1) THEN 
WRITE(10>%)} ° TIME SHOCK CONTACT HE AD 
C TAIL'* 
WRITE{10,"%)} ° ° 
END IF 
WRITE(10,67) T,SIGHMA(1,13,STIGMA(2,1},SIGMA(3,1),SIGMA(4,1} 
67 FORMAT (1X,F12.651X%»F12.651X,F12.6,1X,F12.6,1X,F12.6) 
IF(J.EQ.1) THEN 


WRITE(8,*) ° EXACT VALUES'* 

WRITE(8,%*)} ° TIME SHOCK CONTACT HEAD 
C TAIL’ 

WRITE(8,*} ‘° ' 

END IF 


SKXE =XINIT+VSE*T 
CSXE =XINIT+VCDE XT 
HEWXE =XINIT +VHE ADT 
TEWXE =XINIT+VTAIL*T 
WRITE(8,68) T»SKXE »CSXE »HEWXE » TEWXE 
68 FORMAT (1X%,F12.6,1X,F12.6,1X,F12.6,1X,F12.6,1X,F 12.6) 
RETURN 
END 


SUBROUTINE TIME(N,QQ,RR>S>»,DELT >H} 


FOIE HE HE FETE IE HE IE FE IE HE HE HE IE IE IE IE HEHE HE IE IE IE IE IEIE HE IE IEIE JE IE TE IE TE ETE IE IE HE HE IE IEE HEE 


% % 
% CALCULATE TIME STEP SUBROUTINE * 
¥% % 


Fe HH EIE HE IE ESEE HE HE HEI HE IE IE IE IE IE IE IE IE IE IE IE IEE EE HE IE IE IE EI IEE HE IE IEE HE HIE HM 
-------- NEW VARIABLE DEFINITIONS ------------ 
TMIN - RUNNING VALUE OF THE MINIMUM TIME STEP 


INTEGER N>TI 

DIMENSION QQ(N),RR(N),S(N) 

DOUBLE PRECISION H,A,QQ,RR>S,DELT»TMIN,Q 

TMIN=2.0000 

DO 21 I=1,N 
A=(QQ( I J-RR(I3)7(2.D00%S(T 3) 
Q=( QQ( I I+RRIIII/2.0D00 
DELT=H/{ DABS( DABS(Q)}+A)) 
IF (DELT.LT.TMIN) THEN 
TMIN=DELT 
END IF 

21 CONTINUE 

DELT=0.99D00*TMIN 

RETURN 

END 


SUBROUTINE TRAK(N,SIGMA,H,QQ>RR»S»G»G1»,G2,DELT,I2,X2,W,AR ,DQ>VS> 
#J, LWPRES »QCS,VHEW,VTEW) 


FHKE HE IEE HE HIE IE FEI IE IE HEHE HEHE HE IE IEE IE HEI HEHEHE EH ETE HE HEHE IEE HE HEHE HHH 


pad % 
* DISCONTINUITY TRACKING SUBROUTINE * 
* * 


FE HIE HE IEIE IEE FEFE FE IE IE HE HEE IEE IE HEHE HE IE HE IE IE HE HE IEE HEI IE HEI HE HEI HE IE HEHEHE HH 
woccee----- VARIABLE DEFINITIONS ------------- 


CSDIR - CONTACT SURFACE DIRECTION, 2 TO THE RIGHT,3 TO THE LEFT 
CSRMN - RIEMANN VARIABLE CHANGE ACROSS A CONTACT SURFACE 
DR = THE RATIO OF THE DENSITY ACROSS A 
SHOCK, B/A( RIGHT ),8/A( LEFT ) 
DREMN - DUMMY VARIABLE 
PR = THE RATIO OF THE PRESSURE ACROSS A 
SHOCK, A/B( RIGHT) ,B/A( LEFT) 
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aaa 


MREIMN - THE MEASURED JUMP IN QQ ACROSS THE SHOCK, 
FROM A TO B. 
EREIMN =- THE JUMP IN QQ ACROSS THE SHOCK CALCULATED ANALYTICALLY 
AS A FUNCTION OF W. 
SAP - ENTROPY TO THE LEFT OF THE SHOCK FOR FLOWS RIGHT 
OR ENTROPY TO THE RIGHT OF THE C.S. FOR FLOWS LEFT 
SBP - ENTROPY TO THE RIGHT OF THE C.S. FOR FLOWS RIGHT 
OR ENTROPY TO THE LEFT OF THE SHOCK FOR FLOWS LEFT 
SHKOIR - SHOCK OIRECTION OF TRAVEL>3 TO THE LEFT,e TO THE RIGHT 
TS - TIME FOR SHOCK TO TRAVEL ONE INTERVAL 
X% > DISTANCE FROM LEFT BOUNOARY TO NOOE 


INTEGER N>I>Y»I2,bL>J»SHKOIR»,CSOIR,LWPRES 

DIMENSION SIGMA(4,2),X2(4),RRIN),QQUN),SON),I2(4) 

DOUBLE PRECISION SIGMA >X25X»H»AB»SA+,SB>,AA,QA,QB,QQA ,QQB ,>RRA,RRB> 
RR »QQ,S»TS»N,0Q,AR»,PR>G,»G1»G2,VS,0ELT ,CSRMN> 
Q,MREIMN,OREMN,EREIMN >WW SAIL »SAZ »SAP » SBP 
AW,QWN,VHEW,VTEW, TIME ,QCS 


aan 


+#4444¢44¢ LOCATING THE UPSTREAM NODE +++++¢++++++ 


D0 10 L=1;% 
SIGMA(L,1)=SIGMA(L,2) 
Y=0 
X=0.000 
I=1 
ll IF (.NOT.CY.EQ.0)) GOTO 10 
IF (SIGMA(L,1).LT.X) THEN 
X2(L =X 
IT2(L)=I 
Y=1 
END IF 
X=X+H 
T=I+l 
GOTO 11 
10 CONTINUE 


----- IF SHOCK HAS LEFT AN OPEN BOUNOARY OUT OF THE TUBE THEN SET ---- 
----- SHOCK TO NEUTRAL ------ 


IF (I2(1).GT.N) THEN 
SIGMA(1,1) = 2.000 
SIGMA(1,2) = 3.000 


W = 1.000 

VS = 1.000 
0Q = 0.000 
AR = 0.000 
PR = 1.000 
OR = 1.000 
GO TO 150 


ELSE IF(I2(1).LT.2) THEN 
SIGMA(1,1) =-1.000 
SIGMA(1,2) =-2.000 


W = 1.000 
VS =-1.000 
0Q@ = 0.000 
AR = 0.000 
PR = 1.000 
DOR = 1.000 
GO TO 150 
ENO IF 


$+44+444¢4¢44+4++ DETERMINING SHOCK SPEED +4+4+4444¢4¢ 


IF(CJ.EQ.1).OR.(1T2(01).EQ@.2).OR.(12(1).EQ.N)) THEN 


----- SHKDIR = 3 IS A SHOCK HEADEO LEFT, AND SHKDIR = 2 IS SHOCK----- 
----- TRAVELING RIGHT----- 


4 


o2P 


TFC LWPRES.EQ.3) THEN 


SHKOIR = 3 
IF (J.EQ.1) THEN 
X2(01) = X201) - H 
Batl) = I2t1) -=1 
X2(2) = X202) - H 
I2(2) = T2(2) - 1 
END IF ' 
GO TO 20 
ELSE 
SHKDIR = 2 
ENDO IF 
GO TO 20 
=<--- IF SHOCK AND CONTACT SURFACE ARE NOT WITHIN THE SAME----- 


ELSE IF (I1201).NE.1202))} THEN 
IFC (SHKOTR.EQ.3).AND.(SIGMA(1,1).EQ.(X2(1)-H))) THEN 


X201} = X201)3 = H 
T2(1) = 1241) - 1 
END IF 
GO TO 20 


IF(SIGMA(1,1).GT.SIGMA(2,1)) THEN 


GO TO 111 

ELSE IF(SIGMA(1,1).EQ.(X201)-H)} THEN 
X201) = A201) - H 
LeU yes Letiy =22 
X202) = X202)} - H 
P2t2)e= I2t2) = 1 

ENO IF 

15 RRAS=RR(T201)-1) 


RRB=RR(I2(1)} 
QQA=QQ(T201)-1)3 
QQB=QQ(T2(1)) 
SA=S(T201)-1) 


QA =(QQA+RRA )/2.000 

QB =(QQB+RRB}/72.000 

AA =(QQA-RRA )/(2.D000*SA } 

DQ =(QB-QA)/AA 

W = DSQRT((0Q**2 )*(0.36000 341.000) - (0Q*0.6000} 
OQ =(-1.D00 }*DQ 

GO TO 110 


ELSE IF(SIGMA(1,1).LT.SIGMA(2,1)) THEN 
GO TO 111 

ELSE 

17 RRA=RRII2Z20 13-1) 

RRB=RR(T2(1)) 
QQA=QQ(T2(1)-1) 
QQB=QQ(Tet1)) 
SB=S(IT2¢1)) 
QA =(QQA+RRA )/2.000 
QB =(QQB+RRB }/72.000 
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AB =(QQB-RRB )/(2.0D00*SB ) 
DQ =(QA-QB )/AB 
WH = DSQRT( (DQ#¥2 )¥(0.36900)+2.D00)} + {D0Q#0.6D000) 


GO TO 110 
ENO IF 
C 
C ----- WITH NC SHOCKX/CONTACT SURFACE INTERACTION THE JUMP IN REIMAN----- 
Cc. +-<=- VARIABLES ARE DETERMINED WITHOUT INTERPOLATION OVER THE ----- 
C <<<== INTERVAL. MREIMN IS THE MEASURED JUMP. EREIMN IS THE ANALYTICAL-- 
C s--- VALUE ==<==< 
C 
20 RRA=RR(T2(1)-1) 
RRB=RR(IT201))} 
QQA=QQ(T2(1)-1) 
QQB=QQ(T2(1)) 
SA=S(T201)-1)} 
SB=S(T2(1)} 
21 AB=(QQB-RRB }/(2.000*SB ) 
AA=(QQA-RRA }/(2.DOOXSA ) 
IF(SHKOIR.EQ.3) THEN 
MREIMN = (RRB-RRA}/AA 
DREMN = DABS(MREIMN) 
ELSE 
MREIMN = (QQA-QQB })/AB 
DREMN = MREIMN 
END IF 
Cc 
C ----— ITERATE FOR PROPER VALUE OF W USING THE QUADRATIC FIT OF THE----- 
Cc. ---_-— REIMAN VARIABLE CHANGE WITH W CURVE. NOTE LEFT MOVING SHOCKS---- 
C -----ARE USED IN THESE EQUATIONS SINCE RRB-RRA/AA=-( QQA-QQB/AB )----- 
Cc 
100 WH= (3.0396408D001-( (DREMN+2. 7574000 3/0 . 286337000 } } 
W=5.513294000-DSQRTI WW) 
DQ=2 .000*( W*N-1.D00 J/( W*(G41.D00 }) 
AR=O0SQRT(2.000*(G-1.000 J¥(1.000+((G-1.000 J¥N¥W/2. 000 } }% 

C (G¥G2*¥NXW-1] 000) )/01G41.000 )*W) 
PR=(2.D00*G/(G+1. 000 ) JxW¥N-( (G-1.000)/(G+1.000}} 
DR=((G-1.000 J¥N¥W+2 .000 J/( (G+1.000 )*WN¥W ) 

EREIMN=0Q+( AR-1.D00 )*G2-( AR¥G1/G }*DLOG( PR¥( DRG } } 
IF (DABS(EREIMN-DABS(MREIMN)}.LT.0.10-5) GO TO 110 
DREMN = (DABS(MREIMN) - EREIMN) + OREMN 
GOTO 100 
c 
2---—— SHOCK VELOCITY DEPENDS ON OIRECTION SHOCK IS TRAVELING ----- 
C----- LEFT IS < 0, AND RIGHT IS > 0 <-e-< 
C 


110 IF (J.EQ.1) THEN 
TIME = 0.000 
SAl = (G1/G)*DLOG((2.D00*G i W**2 }-G+1.000 }/(G4+1.D00 }} 
= G1*DLOG( ((G=-1. D000 }¥( We 2 +2 D0 1G41. D000 )*1 We? 1) 
TFISHKOIR.EQ.2) THEN 


SAP = SB - SAl = SA2 
SBP = SAP 
ELSE 
SBP = SA = SAl1 = SA2 
SAP = SBP 
END IF 
103 IF(SHKOIR.EQ.2) THEN 
CSRMN =((DEXP( (SBP=-SA }/G2 ) }*(SA)-(SBP ) }*AR 
ELSE 
CSRMN =( (DEXP( (SAP-SB )/G2 ) )¥(SBJ-CSAP ) JAR 
ENO IF 


IF(SHKDIR.EQ.3) THEN 
MREIMN =( (RRB-RRA }/AA )+CSRMN 
DREMN = DABS(MREIMN } 

ELSE 
MREIMN =( (QQA-QQB }/AB }-CSRMN 
DREMN = MREIMN 

END IF 
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101 WA=(3.03964038001-( ( DREIMN+2. 7574000 )/0. 286337000 ) ) 
W=5.513294D00-DSGRT (WA ) 
0Q=2 .D00*( W*H-1.D00 )/( W4IG+1.000) ) 
AR=DSQRT( 2.0004(G-1.000 i#( 1.000+( (G-1.000 )¥W*W/2 . D000 ) )* 
Cc (G¥G2¥W*N-1.D00 ))/((G+1.D00 )*W) 
PR=(2.D00*G/(G+1.000 ) )*WxW-( (G-1.000 )/(G+1.000)) 
OR=((G-1.D00 )*¥N¥W+2.000 )/((G+1.D00 )*N¥W ) 
ERE IMN=DQ+( AR-1.000 )*G2-( AR¥G1/G )¥DLOG( PR¥( DRG ) ) 
IF (DABSCEREIMN-DABS(MREIMN)).LT.O.10D-5) GO TO 102 ~ 
DREMN = (DABS(MREIMN) = EREIMN) + OREMN 
GOTO 101 
102 SA1 = (G61/G )xDLOG( (2.DO0xG*( Wx*2 )-G+1.D00)/(G+1.D00) ) 

SA2 = G1*¥OLOG( ((G-1. DOO )*(W**2 )4+2)/((G41.D00 )*(WN**2 ) )) 
IF(SHKOTR.EQ.2) THEN 

SAP = 3B - SAl = SA2 
ELSE 

SBP = SA - SAl - SA2 
END IF 
IF (DABS(SAP-SBP}.LT.0.10-5) GO TO 105 
IF(SHKDIR.EQ.2) THEN 

SBP = (SAP-SBP) + SBP 


ELSE 
SAP = (SBP-SAP) + SAP 

ENO IF 

GO TO 103 
ENO IF 

105 IF(SHKOIR.EQ.3) THEN 

OQ = (-1.000)*DQ 

VS = ((RRA+QQA)*0.5000) - (NAA) 
ELSE 

VS=(QQB+RRB )*0. SDOO+W¥AB 
ENO IF 


TS=H/DABS(VS ) 
lll IF (TS.LT.DELT) THEN 
DELT=0.99000*TS 
END IF 
SIGMA( 1,2 )=VS*¥DELT+SIGMA( 1,1) 


4#+++++DETERMINE CONTACT SURFACE SPEED+++++++++4 


IF(J.EQ.1) THEN 
CSOIR = SHKDIR 
END IF 


a CONTACT SURFACE TRAVELING RIGHT ----- 


150 IF(CSDIR.EQ.2) THEN 


IF(J.EQ.1) THEN 
QB =(QQB+RRB)/2.000 
QA = OQ*AB + QB 
CALL SKJUMP(AB,QB,SB8,AR,DQ,VS,G,G1,W,AA,QA)SA ) 
Q = QA 
GO TO 200 
ENO IF 
IF(X2(2).EQ.X2(1)) THEN 
IF(SIGMA(2,1).LE.SIGMA(1,1)) THEN 
Q = (QQ(I2(2)-1) + RRII2(2)-1)) / 2.000 


ELSE 
Q =(QQ(T2(2)) + RR(T2(2))) 7 2.000 
ENO IF 
ELSE 
QA = (QQ(T2(2)-1)} + RRII2(2)-1)) 7 2.000 
QB = (QQ(I2(2)) + RRIT2(2))) 7 2.000 
Q = (QA + QB)/72.000 
ENO IF 


los 
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ELSE IF(CSDIR.EQ.3) THEN 
=--<= CONTACT SURFACE TRAVELING LEFT ----- 


IF(J.EQ.1) THEN 
QA =(QQA+RRA)/2.D00 
QB = DQXAA + QA 
CALL SKJUMP(AA,QA,SA,AR,DQ,VS,G,G1,WN,AB ,QB ,SB) 
Q = QB 
GO TO 200 
END IF 
IF(X202).EQ.X%201)) THEN 
IF(SIGMA(2,1).GE.SIGMA(1,1)) THEN 
Q = (QQ(T2(2)) + RR(T2(2))) 7 2.000 
ELSE 
Q = (QQET2(2)-1) + RROIT2(2)-1)) 7 2.000 
END IF ; 
ELSE IF(SIGMA(2,1).EQ.(X202)-H)) THEN 
R202) = X202) - H 
i2t2) = [2(2) = 1 
Q = (QQ(T2(2))+RRIT202))) 7 2.000 


ELSE 
QA = (QQC(I2(62)-1) + RR(I2(¢2)-1)) 7 2.000 
QB = (Q@Q(I2(2)) # RROII2(2))) 7 2.000 
Q = (QA + QB)/2.000 
END IF 
END IF 
200 SIGMA( 2,2) = DELT * Q + SIGMA( 251) 
QCS=Q 


==--- CALCULATE EXPANSION WAVE SPEED ------ 


TIME = TIME+DELT 
IF(J.EQ.3)} THEN 
AW=(QQ( (N41 0/72 J=-RROON4+1 072 ))702.D00%S((N4+1)72)) 
QH=Q 
IF((CSDIR).EQ.2) THEN 
VHEW=-( AW) 
VTEW= QW-AW 
ELSE 
VHEW= AW 
VTEW= QN+AW 
END IF 
END IF 
TF(J.EQ.3) THEN 
SIGMA(3,2) SIGMA(3,1) + VHEW*TIME 
SIGMA(4,2) SIGMA(G,1) + VTEW*TIME 
ELSE IF((CSDIR.EQ.2).AND.(SIGMA(3,1).GT.0.D00)) THEN 
SIGMA(3,2) SIGMA(3,1) + VHEW*DELT 
SIGMA(4,2)} SIGMA(G,1) + VTEW*DELT 
ELSE IF((CSDIR.EQ.3).AND.(SIGMA(3,1).GT.1.000)) THEN 
SIGMA(3,2) SIGMA(3,1) + VHEWN*DELT 
SIGMA(G,2) SIGMA(G,1) + VTEWXDELT 


END IF 
RETURN 
END 


SUBROUTINE SWEEP(N,H,SIGMA,QQ,RR»S,DELT,EE,Q,A »NENQQ »NEWRR NERS > 
CI2,G2,J,LNODE ,RNODE ,» HALT, LBNORY ,»LBOPRS»LBDPR,LBOTR,LBDDR,QLBD>G, 
CG1,RBNORY »RBOPRS,RBOPR ,>RBOTR»,RBOOR,QRBD >BORY »SK ) 


o.¢.2.2,.¢.0.2.5.5,2, 4.5.2.2, 4,.5,5,5,2,5,5,2,0,2,5,0,¢,5,5, 5,5, 5,.0,5,5,5,. 5,5, 5,2. 0.0.5.5, 5,5,5: 


x x 
ad SPACE SWEEPING SUBROUTINE * 
* * 


3 IE FE FETE IE IEE FE HE IEE IE IE IE IE IE IEE IE IE IEE IE IEE IE IEE IE IE IE IEE IEE IE IE IE IE IE IEE EE EE 


------------ VARIABLE DEFINITIONS ------------ 


MOaAaANNANAANNANANNNNNANNAaNNAANAANAHNgANAAANAAaNAAnANn 


AAVG - AVERAGE SPEED OF SOUND 
CNTACT = 3 DIGIT VARIABLE DENOTING CONTACT SURFACE 

LOCATION, CIRECTION OF TRAVEL, AND IF IT CROSSES A NODE 
DELQQH - CHANGE IN GQ FROM I TO I+1 
DELQ@QL - CHANGE IN @Q FROM I-1 TO TI 
DELRRH = CHANGE IN RR FROM I TO I+1 
DELRRL - CHANGE IN RR FROM I-1 TO TI 
DELSH - CHANGE IN S$ FROM I TO I+1 
DELSL = CHANGE IN S$ FROM I-1 TOT 
OELAH = CHANGE IN A FROM I TO I+41 
DELAL - CHANGE IN A FROM I-1 TO T 
DELQH = CHANGE IN Q FROM I TP I+1 
DELQL = CHANGE IN @ FROM I-1 TO I 
DELX - INTERPOLATION DISTANCE (LMDXDELT ) 
DOLTA*% = PREFIX WHICH INDICATES THE SPATIAL 

CHANGE IN *% FOR ONE TIME STEP. 
INTEG(K )=- RESULT OF INTEGRATING 2(K) 
¥*INT - VALUE OF *# INTERPOLATED BETWEEN NODES 

ON THE CURRENT TIME LEVEL. 
LXX = NODE DEFINING THE LEFT INTERVAL 
*PRIM(K J= SUFFIX WHICH INDICATES THE SPATIAL 

OERIVITIVE OF * AT THE CURRENT TIME LEVEL. 
RXX = NODE DEFINING THE RIGHT INTERVAL 
SAVG - AVERAGE ENTROPY 
SHOCK - 3 DIGIT VARIABLE OENOTING SHOCK 

LOCATION,DIRECTION OF TRAVEL, ANDO IF IT CROSSES A NODE 
*#*HSTEP =- THE CHANGE IN TIME OF ** AT A NODE 

USEO TO STEP UP TO THE NEXT TIME LEVEL 
x - LOCATION IN SPACTAL PLANE (I-1)#H 
Z(K) - RIGHT SIOE OF THE K'TH EQUATION. 


INTEGER I,RXX,LXX,SHOCK ,CNTACT ,I2,J,LNODE,RNODE,HALT, 
N»LBNORY ,LBOPRS ,RBNORY ,RBOPRS »SK,BORY 
DIMENSION SIGMA(4,2),S(N),Q(N),40N),I2(¢),QINT(3 ),AINT(3),2(3), 
NEWQQ(N),NEWRR(N),NEWNS(N),INTEG(3),APRIM(3),QPRIM(3), 
LNODE(4),RNODE(4),AAVG(3)},RRIN),QQ(UN) 
DOUBLE PRECISION AAVG,SAVG,G2 5X,H,SIGMA,QQ,RR»S,Q,A,G6>,G1,5 
DELQQH ,DELQQL ,DELRRH,DELRRL,DELSH,DELSL,> 
DELAH,OELAL ,DELQH,DELQL,DELT, 
QINT ,AINT ,QQINT ,RRINT ,SINT,EE> 
NEWQGQ ,NEWRR ,NEWS ,»LBDPR,LBOTR,LBOOR,QLBO, 
RRSTEP,SSTEP ,INTEG,Z,OLTAQQ,QQSTEP, 
OLTARR,OLTAS ,APRIM,QPRIM, 
RBOPR,RBOTR,RBOOR,QRBD 
COMMON AR,0Q,VS5,W 


ono no 


aangaaan 


OS COMPUTE VELOCITY AND SPEED OF SOUND AT EACH NODE ----- 


0O 10 I= 1,N 
Q(I) = (QQ(I) + RR(II)) 7 2.0000 
= ((QQ(T) - RRII)) 702.0000 * S(T))) 


=~ a 


10 CONTINUE 


sence ADVANCE LEFT BOUNDARY TO NEW TIME STEP ----- 


BORY = 3 
IF(I2(1).EQ.2) THEN 
SK = 1 
ELSE IF(SIGMA(1,2).EQ.-2.000) THEN 
SK = 3 
ELSE 
SK = 0 
ENO IF 
CALL BONORY(Q(1),Q(2),QLB0,A(1),A02),QQ(1),QQ(2),RR(1),RRI2)>5 
Cc $(1),S(02),H,EE,OELT,LBNORY,LBOPRS,LBOPR,LBDDR,LBOTR, 
Cc J ,NENQQ(1),NENRRO1),NENS(1),6,61,G2,HALT,8DRY »SK } 
Sas AT EACH NODE FROM 2 TO N-l DETERMINE THE BEST ALGORITHM TO USE-- 
———= TO ADVANCE THAT NODE TO THE NEXT TIME STEP ----- 
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I=2 
11 IF(I.EQ@.N) GO TO 1200 
X=FLOAT( I-12 )*H 
DELQQH = QQ(T+1) = QQ(T) 


DELQQL = QQ(T) - QQ(I-1) 
DELRRH = RR(I+1) = RRITI) 
DELRRL = RR(II) - RR(II-1) 
DELSH = S(I+#1) = S(T) 
DELSL = S(I) - S(I-1) 
DELAH = A(T+1) = A(T) 
DELAL = A(T) - AlTI-1) 
DELQH = Q(I+#1) - Q(T) 
DELQL = Q(I) - Q(I-1) 


IF (1201).EQ.RXK) THEN 


SHOCK = 20D 
GO TO 20 
ELSE IF (I12(1).EQ.LXX) THEN 
SHOCK = 30D 
GO TO 2D 
EUSE 
SHOCK = 1DD 
END IF 
GO TO 3D 


2D IF (SIGMA(1,1).LT.SIGMA(1,2)) THEN 
SHOCK = SHOCK + 20 
EUSE 
SHOCK = SHOCK + 3D 
END IF 


7-oce DETERMINE IF SHOCK CROSSES A NODE IN THIS TIME INTERVAL ----- 


IF (SHOCK.EQ.220) THEN 
- IF (SIGMA(1,2).GE.(X+#H)) THEN 
SHOCK = SHOCK + 1 
ELSE 
SHOCK = SHOCK + 2 
ENO IF 
ELSE IF (SHOCK.EQ.230) THEN 
IF (SIGMA(1,2).LE.X) THEN 
SHOCK = SHOCK + 1 
ELSE 
SHOCK = SHOCK + 2 
END IF 
ELSE IF (SHOCK.EQ.320) THEN 
IF (SIGMA(1,2).GE.%) THEN 
SHOCK = SHOCK + 1 
ELSE 
SHOCK = SHOCK + 2 
END IF 
ELSE IF (SHOCK.EQ.330) THEN 
IF (SIGMA(1,2).LE.(X-H)) THEN 
SHOCK = SHOCK + 1 
ELSE 
SHOCK = SHOCK + 2 
END IF 
END IF 
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——--— TEST FOR CONTACT SURFACE ----- 


30 IF (12(2).EQ.RXX) THEN 

CNTACT = 200 
GO TO 40 

ELSE IF (12(2).EQ.LXX)} THEN 
CNTACT = 300 
GO TO 40 

ELSE 
CNTACT = 100 

END IF 

GO TO 50 


GO IF (SIGMA(2,1).LT.SIGMA(2,2)}) THEN 
‘CNTACT = CNTACT + 20 
ELSE 
CNTACT 3 CNTACT + 30 
END IF 
c 
C ----- DETERMINE IF CONTACT SURFACE CROSSES A NODE DURING THIS TIME --- 
C ----- INTERVAL ----- 


IF (CNTACT.EQ.220) THEN 
IF (SIGMA(2,2).GE.(X+H)) THEN 
CNTACT = CNTACT + 1 
ELSE 
CNTACT = CNTACT + 2 
ENO IF 
ELSE IF (CNTACT.EQ.230) THEN 
IF (SIGMA(2,2).LE.4) THEN 
CNTACT = CNTACT + 1 
ELSE 
CNTACT = CNTACT + 2 
END IF 
ELSE IF (CNTACT.EQ.320) THEN 
IF (SIGMA(2,2).GE.X%) THEN 
CNTACT = CNTACT + 1 
ELSE 
CNTACT + 2 


CNTACT 
END IF 
ELSE IF (CNTACT.EQ.330) THEN 
IF (SIGMA(2,2).LE.(X-H))} THEN 
CNTACT = CNTACT + 1 
ELSE 
CNTACT = CNTACT + 2 
END IF 
ENO IF 


@--CHECK IF EITHER A SHOCK OR CONTACT SURFACE WITHIN H OF THIS NODE--- 
@--DETERMINE PROPER ALGORITHM TO USE FOR CALCULATING EXTENDED REIMAN-- 
---VARIABLE CHANGE ALONG CHARACTERISTICS AT THIS NODE----- 


aanNan 


50 IF (SHOCK.EQ@.100.OR.CNTACT.EQ.100) THEN 
“--NEITHER A SHOCK NOR A CONTACT SURFACE EXIST NEAR NODE--- 
IF (SHOCK.EQ.100.AND.CNTACT.EQ.100) THEN 


---TEST FOR SUBSONIC OR SUPERSONIC FLOW--- 


ao0o”0n aan 


IF (DABS(Q(I)).LT.A(I})) THEN 
IF(T.EQ.(12(1)-2)) THEN 
IF(J.EQ.1) THEN 
IF(SIGMA(1,1).LE.SIGMA(1,2)) THEN 
GO TO 200 
; ENO IF 
END IF 
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END IF 
IF(I.EQ.¢I12¢1)41)) THEN 
IF(J.EQ.1) THEN 
IF(SIGHA(1,1).GE.SIGMA(1,2)) THEN 


GO TO 300 
END IF 
END IF 
END IF 
GO TO 100 
ELSE 
IF (SIGMA(2,1).LE.SIGMA(2,2)) THEN 
GO TO 200 
ELSE 
GO TO 300 
END IF 
END IF 


END IF 
---SHOCK OR CONTACT SURFACE ON LEFT, HEADED RIGHT, NO NODES CROSSED-- 


IF (SHOCK .EQ.322.0R.CNTACT.EQ.322)} THEN 
IF (DABS(Q(I)).LT.A(I}) THEN 
GO TO 300 
ELSE 
GO TO 500 
END IF 
END IF 


---SHOCK OR CONTACT SURFACE ON LEFT, HEADED LEFT, NO NODES CROSSED-- 


IF (SHOCK .EQ.332.OR.CNTACT.EQ.332) THEN 
GO TO 300 
END IF 


~--SHOCK OR CONTACT SURFACE ON LEFT, HEADED RIGHT, NODE IS CROSSED--- 


IF (SHOCK.EQ.321.0OR.CNTACT.EQ.321) THEN 
GO TO 400 
END IF 


---SHOCK OR CONTACT SURFACE ON LEFT, HEADED LEFT, NODE IS CROSSED--- 


IF (SHOCK.EQ.331.OR.CNTACT.EQ.331) THEN 
GO TO 300 
END IF 


“--SHOCK OR CONTACT SURFACE ON RIGHT, HEADED RIGHT, NO NODES CROSSED-- 


IF (SHOCK.EQ.222.0R.CNTACT.EQ.222)} THEN 
GO TO 200 
END IF 


“--SHOCK OR CONTACT SURFACE ON RIGHT, HEADED LEFT, NO NODES CROSSED-- 


IF (SHOCK.EQ.232.0R.CNTACT.EQ.232) THEN 
IF (DABS(Q(I)).LT.ACI)) THEN 
GO TO 200 
ELSE 
GO TOC 500 
END IF 
END IF 


---SHOCK OR CONTACT SURFACE ON RIGHT, HEADED RIGHT, NODE CROSSED--- 
IF (SHOCK .EQ.221.OR.CNTACT.EQ.221) THEN 
GO TO 200 
END iF 


---SHOCK OR CONTACT SURFACE ON RIGHT, HEADED LEFT, JUMPS NODE--- 
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GO TO 400 
END IF 
Cc 
C ---BRANCH HERE IF SHOCK TO RIGHT OF CONTACT SURFACE ----- 
C ---DETERMINE PROPER ALGORITHM TO USE FOR CALCULATING EXTENDED REIMAN-- 
C ---VARIABLE CHANGE ALONG CHARACTERISTICS AT THIS NODE----- 
c 
IF (SIGMA(1,1).GT.SIGMA(2,1)) THEN 
Cc 
C ---SHOCK ON LEFT, HEADED RIGHT, NO NODE JUMPED--- 
Cc 
IF (SHOCK.EQ.322) THEN 
IF (CNTACT.EQ.322) THEN 
IF (DABS(Q(I)).LT.A(T)) THEN 
GO TO 300 
ELSE 
GO TO 500 
END IF 
ELSE IF (CNTACT.EQ.332) THEN 
GO TO 800 
ELSE IF (CNTACT.EQ.331) THEN 
GO TO 800 
ELSE 
GO TO 900 
END IF 
END IF 
Cc 
C =---SHOCK ON LEFT, HEADED LEFT, NO NODE JUMPED--- 
C 
IF (SHOCK.EQ.332) THEN 
IF (CNTACT.EQ.322) THEN 
GO TO 300 
ELSE IF (CNTACT.EQ.332) THEN 
GO TO 600 
ELSE IF (CNTACT.EQ.331) THEN 
GO TO 600 
ELSE IF (CNTACT.EQ.321) THEN 
GO TO 700 
ELSE 
GO TO 900 
ENO IF 
END IF 
Cc 
C s--SHOCK ON LEFT, HEADED RIGHT, JUMPS NODE--- 
C 
IF (SHOCK.EQ.321) THEN 
IF (CNTACT.EQ.322.OR.CNTACT .EQ.321) THEN 
GO TO 400 
ELSE IF (CNTACT.EQ.332.OR.CNTACT.EQ.331) THEN 
GO TO 800 
ELSE 
GO TO 900 
END IF 
END IF 
C 
C ---SHOCK ON LEFT, HEADED LEFT, JUMPS NODE--- 
C 
IF (SHOCK.EQ.331), THEN 
IF (CNTACT.EQ.322.OR.CNTACT.EQ.321) THEN 
GO TO 700 
ELSE iF (CNTACT.EQ.332.OR.CNTACT.EQ.331) THEN 
GO TO 600 
ELSE 
GO TO 900 
END IF 
END IF 
Cc ’ 


C ---SHOCK ON RIGHT, HEADED RIGHT, NO NODE JUMPED--- 
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Cc 


IF (SHOCK.EQ.c22} THEN 


END IF 


IF (CNTACT.EQ.222) THEN 
GO TO 200 

ELSE IF (CNTACT.EQ.232.OR.CNTACT.EQ.231) THEN 
GO TO 800 

ELSE IF (CNTACT.EQ.321.OR.CNTACT.EQ.322) THEN 
GO TO 400 

ELSE IF (CNTACT.EQ.332.OR.CNTACT.EQ.331) THEN 
GO TO 800 

EESE 
GO TO 900 

ENO IF 


C ---SHOCK ON RIGHT, HEAOQEO LEFT, NO NOOE JUMPEO--- 


Cc 


Cc 


IF (SHOCK.EQ.232) THEN 


ENO IF 


IF (CNTACT.EQ.222) THEN 
IF (SIGMA(1,2).LT.SIGMA(2,2)) THEN 


GO TO 700 
ELSE 
GO TO 200 
ENO IF 
ELSE IF (CNTACT.EQ.231.OR.CNTACT.EQ.232) THEN 
GO TO 600 
ELSE IF (CNTACT.EQ.221) THEN 
GO TO 700 
ELSE IF (CNTACT.EQ.322) THEN 
GO TO 400 
ELSE IF (CNTACT.EQ.332.0R.CNTACT.EQ.331) THEN 
GO TO 600 
ELSE IF (SIGMA(1,2}.LT.SIGMA(2,2)} THEN 
GO TO 710 
ELSE 
GO TO 400 
ENO IF 


C ---SHOCK ON RIGHT, HEAQEO RIGHT, JUMPS NOOE--- 


Cc 


Cc 


IF (SHOCK .EQ.221) THEN 


END IF 


IF (CNTACT.EQ.221.0R.CNTACT.EQ.222) THEN 
GO TO 200 

ELSE IF (CNTACT.EQ.231.OR.CNTACT.EQ.232) THEN 
GO TO 800 

ELSE IF (CNTACT.EQ.321.OR.CNTACT.EQ.322) THEN 
GO TO 400 

ELSE == 
GO TO 800 

ENO IF 


C ---SHOCK ON RIGHT, HEAQEO LEFT, JUMPS NOOE--- 


Cc 


IF (CNTACT.EQ.221.OR.CNTACT.EQ.222) THEN 


GO TO 700 


ELSE IF (CNTACT.EQ.231.OR.CNTACT.EQ.232) THEN 


GO TO 600 


ELSE IF (CNTACT.EQ.322} THEN 


IF (SIGMA(1,2).LT.SIGMA(2,2)) THEN 
GO TO 710 

ELSE 
GO TO 400 

END IF 


ELSE IF (CNTACT.EQ.321) THEN 


GO TO 710 


ELSE IF (OABS(Q(I)).LT.A(I)) THEN ’ 


GO TO 600 
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ELSE 
GO TO 800 
END IF 


---BRANCH HERE IF SHOCK IS TO LEFT OF CONTACT SURFACE---- 
---DETERMINE PROPER ALGORITHM TO USE FOR CALCULATING EXTENDED REIMAN-- 
@“--VARIASLE CHANGE ALONG CHARACTERISTICS AT THIS NODE----- 

ELSE IF (SIGMA(1,1).LT.SIGMA(2,1)) THEN 


~--SHOCK ON LEFT, HEADED RIGHT, NO NODE CROSSED--- 


aan nAanaaAn 


IF (SHOCK.EQ.322) THEN 
IF(CNTACT .EQ.322.0R.CNTACT.EQ.321) THEN 
GO TO 600 
ELSE IF (CNTACT.EQ.332) THEN 
IF (SIGMA(1,2).GT.SIGMA(2,2)) THEN 
GO TO 700 
ELSE 
GO TO 300 
END IF 
ELSE IF (CNTACT.EQ.331) THEN 
GO TO 700 
ELSE IF (CNTACT.EQ.222.0R.CNTACT.EQ.221) THEN 
GO TO 600 
ELSE IF (CNTACT.EQ.232) THEN 
GO TO 400 
ELSE IF (SIGMA(1,2}.GT.SIGMA(2,2)) THEN 
GO TO 710 
ELSE 
GO TO 400 
ENO IF 
END IF 
c 
C ---SHOCK ON LEFT, HEADED LEFT, NO NODE CROSSED--- 
c 
IF (SHOCK.EQ.332) THEN 
IF (CNTACT.EQ.322) THEN 
GO TO 800 
ELSE IF (CNTACT.EQ.332) THEN 
GO TO 300 
ELSE IF (CNTACT.EQ.321) THEN 
GO TO 800 
ELSE IF (CNTACT.EQ.222.OR.CNTACT.EQ.221) THEN 
GO TO 3800 
ELSE IF (CNTACT.EQ.231.OR.CNTACT.EQ.232) THEN 
GO TO 400 
ELSE . 
GO TO 900 
ENO IF 
END IF 
c 
C ---SHOCK ON LEFT, HEADED RIGHT, JUMPS NODE--- 
© 
IF (SHOCK.EQ.321) THEN 
IF (CNTACT.EQ.322.0R.CNTACT.EQ.321) THEN 
GO TO 600 
ELSE IF (CNTACT.EQ.332.0R.CNTACT.EQ.331) THEN 
GO TO 710 
ELSE IF (CNTACT.EQ.222.0R.CNTACT.EQ.221) THEN 
GO TO 600 
ELSE IF (CNTACT.£Q.231) THEN 
GO TO 710 
ELSE IF (SIGMA(1,2).GT.SIGMA(2,2)) THEN 
GO TO 710 
ELSE 
GO TO 400 
END IF 
END IF 
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C 
C ---SHOCK ON LEFT, KEADEO LEFT, JUMPS NODE-~- 
Cc 
IF (SHOCK.EQ.331) THEN 
IF (CNTACT.EQ.322.0R.CNTACT.EQ.321) THEN 
GO TO 800 
ELSE ITF (CNTACT.EQ.331.OR.CNTACT.EQ.332) THEN 
GO TO 300 
ELSE IF (CNTACT.EQ.221.OR.CNTACT .EQ.222) THEN 
GO TO 800 
ELSE 
GO TO 400 
END IF 
END IF 
Cc 
C ---SHOCK ON RIGHT, HEADED RIGHT, NO NODE CROSSED--- 
Cc 
IF (SHOCK.EQ.222) THEN 
IF (CNTACT.EQ.222.0R.CNTACT.EQ.221) THEN 
GO TO 600 
ELSE IF (CNTACT.EQ.231) THEN 
GO TO 710 
ELSE IF (CNTACT.EQ.232) THEN 
IF (SIGMA(1,2).GT.SIGMA(2,2)) THEN 


GO TO 700 
ELSE 
GO TO 200 
ENO IF 
ELSE 
GO TO 900 
ENO IF 


END IF 
Cc 
C ---SHOCK ON RIGHT, HEADED LEFT, NO NODES CROSSED--- 
Cc 
IF (SHOCK.EQ.232) THEN 
IF (CNTACT.EQ.222.OR.CNTACT .EQ.221) THEN 
GO TO 800 
ELSE IF (CNTACT.EQ.232) THEN 
IF (DABS(Q(I}3.LT.ACT 33 THEN 
GO TO 200 
ELSE 
GO TO 500 
END IF 
ELSE 
GO TO 900 
END IF 
END IF 
Cc 
C ---SHOCK ON RIGHT, HEADED RIGHT, JUMPS NODE--- 
Cc 
IF (SHOCK.EQ.221) THEN 
IF (CNTACT.EQ.221.OR.CNTACT.EQ.222) THEN 
GO TO 600 
ELSE IF (CNTACT.EQ.231) THEN 
GO TO 710 
ELSE IF (CNTACT.EQ.232) THEN 
GO TO 700 
ELSE 
GO TO 900 
END IF 
END IF 
Cc 
C ---SHOCK ON RIGHT, HEADEO LEFT, JUMPS NODE--- 
Cc 
IF (CNTACT.EQ.221.OR.CNTACT.EQ.222) THEN 
GO TO 800 
ELSE IF (CNTACT.EQ.231.OR.CNTACT .EQ.232) THEN 
2 GO TO 400 | 
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c 


ELSE 
GO TO 900 
END IF 


C -=-SHOCK AND CONTACT SURFACE ARE AT THE SAME LOCATION AFTER TIME--- 
C ---ZERO--- LG 


c 


aAanaAn 


ELSE IF((SHOCK.EQ.222).AND.(CNTACT.EQ. 222) JTHEN 


SHOCK = 321 
CNTACT = 321 
GO TO 400 


ELSE IF( (SHOCK .EQ.322).AND.(CNTACT .EQ.322) )THEN 
IF(DABS(Q(I)).LT.ACI)) THEN ; 
GO TO 300 
ELSE 
GO TO 500 
END IF 
ELSE IF( (SHOCK .EQ.332).AND.(CNTACT.EQ.332)) THEN 
SHOCK = 231 
CNTACT = 231 
GO TO 400 
ELSE IF( (SHOCK .EQ.232).AND.(CNTACT.EQ.232)) THEN 
IF(DABS(Q(I)).LT.ACI)) THEN 


GO TO 200 
ELSE 
GO TO 500 
END IF 
ELSE 
GO TO 720 
END IF 


~--CALL CONDITION SUBROUTINE WHICH CONTAINS ALGORITHM THAT IS--- 
m--NUMERICALLY 8EST SUITED FOR THE SITUATION AT NODE I--- 


100 CALL COND1(Q(TI),Q(0I+1),A(T),A0CT*1),RRII),QQ(IT),S(I),DELQQL, 


Cc 
Cc 


CELRRH,DELSH,DELSL,DELQH,DELAH,OELQL ,DELAL,H,EE, 
DELT »,QQINT ,RRINT ,SINT ,QPRIM,APRIM,AINT ,Q(I-1),A(I-1)) 


GO TO 1000 


200 CALL COND2(Q(1I),Q(IT-1),A(I),A(T-1),RRIT),QQ(I),S(I),DELQQL, 


Cc 
Cc 


OELRRL,DELSL,DELQL,DELAL,DELT,H,EE,QQINT ,RRINT >SINT > 
QPRIM,APRIM,AINT ) 


GO TO 1000 


300 CALL COND3(Q(I),Q(I+1),A(T),A0T+1),RRIT),QQ(1T),S(1),D0ELQQH, 


Cc 
Cc 


DELRRH,DELSH ,DELQH,DELAH,DELT,H,EE,QQGINT,RRINT,SINT, 
QPRIM,APRIM,AINT ) 


GO TO 1000 


400 CALL COND4(T,SHOCK ,CNTACT,J,LNOOE ,RNODE ,QQSTEP,RRSTEP,SSTEP ) 


GO TO 1100 


500 CALL CONDS( Q(T), Q(IT-1),Q(I4+1),A0T),A(T-1),A(CT4+1),RRIT),QQIT), 


Cc 
Cc 
C 


GO TO 1000 


S(I},DELQQL ,DELQQH,DELRRL ,DELRRH,DELSL ,DELSH,OELQL, 
DELQH,DELAL ,DELAH,H,EE ,DELT,QQINT »RRINT »,SINT >AINT, 
QPRIM,APRIM ,SHOCK ,CNTACT ) 


600 CALL COND6( SHOCK ,CNTACT »HALT ) 


QQSTEP=0.00 
RRSTEP=0.D0 
SSTEP=0.00 
GO TO 1100 


700 CALL CONO7(SHOCK ,CNTACT ,DELT,SIGMA,I2,I,H,HALT ,QUT )) 


QQSTEP=0.00 
RRSTEP=0.D00 
SSTEP=0'.00 
GO TO 1100 
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710 CALL COND7N( SHOCK ,CNTACT ,DELT >SIGMA,I2,I>H» HALT ,Q(I),X3 


QQSTEP=0.D0 

RRSTEP=0.D0 
SSTEP=0.D0 

GO TO 1100 


720 CALL COND7S(SIGMA ,HALT ,»SHOCK ,CNTACT ) 
QQSTEP=0.00 
RRSTEP=0.D0 
SSTEP=0.D0 
GO TO 1100 


800 CALL CONOD8( SHOCK ,CNTACT >HALT ) 
QQSTEP=0.D0 
RRSTEP=0.00 
SSTEP=0.00 
GO TO 1100 


900 PRINT * »'AN IMPOSSIBLE SITUATION EXISTS - ERROR‘ 
QQSTEP=0.D0 
RRSTEP=0.D0 
SSTEP=0.D0 
HALT = 1 
GOTO 1100 


wwoccer- CALCULATE DLTA QQ, DLTA RR & DLTA S 


1000 DLTAQQ=QQINT-QQ(T) 


DLTARR=RRINT-RRIT } 
OLTAS=SINT-S(T ) 


~----- CALCULATE Z(K)'S ----------- 


AAVG(1 )=( AINT(1 J4A0T 172.0000 
AAVG(3 J=(AINT(3 J4#A(T } 372.0000 
AAVG( 2 }=0.0000 

SAVG = (SINT#S(1I))72.000 


Z2(1)=-(1.0000/G2 )*AAVG( 1 ¥CSAVG-G2 }*(QPRIM( 1 }-G2*APRIM(1)) 
Z(3 }=(1.0D00/G2 }¥AAVG( 3 }*( SAVG-G2 }*( QPRIM( 3 }+G2¥APRIM( 3} } 


2{2)=0.0000 


INTEG( 2 )=0.0000 
INTEG(1 )=Z(1 DXDELT 
INTEG( 3 }=2(3 )¥DELT 


QQSTEP=OLTAQQtINTEG( 1) 
RRSTEP=DLTARR+INTEG( 3 ) 
SSTEP=DLTAS+tINTEG( 2 } 


1100 NEWQQ(T }=QQ(T )+QQSTEP 


NEWRR( TI J=RRCI J+RRSTEP 
NEWS(T J=SCI }+SSTEP 


--------- GO TO NEXT NODE ------------- 


T=I+1 
GOTO 11 


1200 CONTINUE 


BORY = 2 
IF(I2(01).EQ.N)} THEN 
SK = 1 
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1210 IF 


ELSE IF(SIGMA{1,23.EQ@.3.000) THEN 
SK = 3 

ELSE 
SK = 0 

END IF 

CALL BONDRY( Q(N),Q(N-1),QRBO,A(N),A0N-1),QQ(N),QQ(UN-1),RRIN), 

RR(N-1),S(N),SCN-1),H,EE ,DELT »RBNDRY ,RBOPRS»RBOPR,RBOOR >» 

RBDTR »J,NENQQ(N) »,NEWRR(N ),NEWS(N),G,G1,G2,HALT,BORY»SK } 


------ UPDATE THE VARIABLES ---------- 


(I.EQ.N4#1) GOTO 1220 
RR(I J=NEWRROT ) 
QQUT J=NENQQ(T ) 
SUI J=NEWS(T ) 


1220 CONTINUE 


T=I+1 

GOTO 1210 

RETURN 

END 

3 EE IEE SESE FE IE IEE SIE IE IESE BE SEBEL HE BE 2G C3 BE MIE IE SEE AEE IEE BE IE IEE AE IES MEM CIEE LE EHC 3EIE 
% % 
% CONDITION SUBROUTINES 1-8 % 
% % 


HEIE IE HE FETE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE HE IE IE IE IE IE IE IE JE EIE IE IE IE IE IE IE IE HE IE IEIE IE IE IE IE IE IE IE IE IEIE IE HEHE 


SUBROUTINE COND1(QI ,QIP1,AI,AIP1»RR»QQ,S,DELQQL ,DELRRH,DELSH,» 
DELSL,DELQH,DELAH,DELQL ,DELAL ,H,EE,DELT,QQINT, 
RRINT »SINT ,»,QPRIM,APRIM,AINT »QIM1,AIM1 ) 


JEIE 3G FETE HEI FE IE 3G IE IE IE IE IEE IE ETE IE FEE IE IE IE IE IE IE IE IE IE IE IE HE IE IEE FETE IE 3G IE 3 IE IE IE IE 36 3E IE FETE 3E IE IE HE 


% % 
ad SUBROUTINE CONDITION 1 * 
* % 


HEHEHE IE HE IE IE IE HE IEE HE IE IE IE IE IE IE IE IE IEE OIETE 3G IE IE IE HE HE IE IE IE IE IE IEIE HEHE IE IE IE IE IE IE FE IE IE IE IE IE IE IE IEE 


———— USED WHEN NO CONTACT SURFACES NOR SHOCKS WITHIN H OF NODE----- 
a= AND FLOW IS SUBSONIC----- 


----- c 


ALCULATES QQINT,RRINT SINT ,QPRIM, APRIM----- 


=<22- USES FORWARD AND BACKWARD DIFFERENCE SCHEMES----=- 


AI - A(T). 

AIM1 - A(I-2) 

AIP] _ A(I+1) 

E(K) - ACTUAL ERROR IN CHARACTERISTIC SLOPE CALCULATION. 
LMD § - SLOPE OF THE CHARACTERISTICS (Q+A,Q,Q-A) 

QI - Q(T) 

QIM1 - Q(I-1) 

QIP1 - Q(I+1) 


aAaM 


DIMENSION LMD(3),DELX(3),QINT(3),AINT(3),E(3},QPRIM(3),APRIM(3 ) 
INTEGER K 
DOUBLE PRECISION QI,QIP1,AI ,AIP1,RR»,QQ,S,AIM1 ,QIM1, 
DELQQL ,DELRRH,DELSH,DELSL,DELQH,DELAH,DELQL, 
DELAL ,H,EE,DELT,LMD,DELX,QINT ,AINT>E, 
QQINT »RRINT »>SINT »,QPRIM,APRIM 


INITIAL ESTIMATE OF CHARACTERISTIC SLOPES----- 


LMD(1) = QIM1 + AIM1 
LMD(2) = QI 
LMD(3) = QIP1 - AIP1 
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10K =1 


20 IF (K.LT.4) TREN 
CDELX(K) = DELT * LMD(K) 
IF (LMO(K).LT.0.0000) THEN 


QINT(K) = QI - (DELX(K) * CELQH / H) 
AINT(K) = AI + (DELX(K) * DELAH 7/7 H) 
ELSE 
QINT(K) = QI - (DELX(K) ® DELQL / H) 
AINT(K) = AI - (DELX(K) * DELAL / H) 
END IF 
K=K #1 
GO TO 20 


END IF 


==s— = CALCULATE ERROR BETWEEN ESTIMATED SLOPE AND NEW SLOPE----- 
Salad FROM NEW INTERPOLATED VALUES----- 


E(t?) 
E(2) 
E(3) 


LMD(1) 
LMD( 2) 
LMD(3) 


DABS(LMD(1) - (QINT(1) # AINT(1))) 
DABS(LMD(2) - (QINT(2))) 
DABS(LMD(3) - (QINT(3) = AINT(3))) 


QINT(1) # AINT(1) 
QINT(2) 
QINT(3) - AINT(3) 


=<=-- COMPARE ERROR TO ERROR TOLERANCE LEVEL, ITERATE TIL MET----- 


IF ((E(1).GT.EE).OR.(E(2).GT.EE).OR.(E(3).GT.EE)) GO TO 10 


socce CALCULATE LINEARLY INTERPOLATED VALUES OF REIMAN----- 
a VARIABLES AND MODIFIED ENTROPY AT POINT A----- 


QQINT = QQ - (DELX(1) * (DELQQL / H)) 
IF (LMD(2).LE.0.0000} THEN 
SINT = S$ - (DELX(2) * (DELSH 7 H)} 


ELSE 


SINT = S = (DELX(2) * (DELSL 7 H)) 


ENO IF 


RRINT = RR - (DELX(3) * (DELRRH 7 H)) 


“-~ CALCULATE SPATIAL DERIVATIVES --- 


QPRIM(13 
QPRIM( 2) 
QPRIM(3 } 
APRIM(1) 
APRIM( 2) 
APRIM(3) 
RETURN 
END 


Te | | | | 


(QI-QIM13/H 
0.000 
(QIP1-QI )/H 
(AI-AIM1)/H 
0.000 
(AIP1-AI )/H 


SUBROUTINE COND2(QI ,QIM1,AI ,AIM1,RR,QQ,S,DELQQL,DELRRL,DELSL, 


DELQGL ,DELAL,DELT,H,EE ,QQINT ,RRINT >SINT ,QPRIM,> 
APRIM,AINT ) 


FE HE IE IE HE IESE FEE HE HE HE IE IE IE IE HE IE IE IE IE IE IEE HE HE IE IEE FE IE FEE HE IE IE IE IE IE HE IE TE IEE FE IE HE IE DE FE IE EE 3 3G 3 


% 
* 
% 


x 
SUBROUTINE CONDITION 2 * 
% 


FETE IEE HE IEE IE IE IE HE IE IE IE IE IE IE IE IE IE IE IE IE IE ETE HE IE FE IE IE FE IE IE IE HE IE IE IE IE IE IE IE IE IE IE IE IE IE IE HE IE IE IE HE HE 


--~ BACKWARD DIFFERENCE ALGORITHM FOR CALCULATING --- 
--- RRINT,QQINT,SINT,APRIM,QPRIM,AINT --- 


OTMENSION LND(3),DELK(3},QINT(3),AINT(3),E(3},QPRIM(3),APRIM( 3} 


INTEGER K 


DOUBLE PRECISION QI ,QIM1,AI,AIM1,RR,QQ,S,DELQQL,DELRRL,DELSL, 


DELQL,DELAL,DELT,H,EE,AINT ,QINT,LMD,DELX>E, 
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Cc QQINT »>RRINT ,>SINT ,QPRIM,APRIM 


wo----INITIAL ESTIMATE GF CHARACTERISTIC SLOPES----- 


LMD(1) = QIM1 + AIML 
LMD(2) = GI 
LMD(3) = QI - AT 


w=---CALCULATE LINEARLY INTERPOLATED VALUES OF @ AND A----- 


10K =1 

20 IF (K.LT.4) THEN 
DELX(K) = DELT * LMD(K) 
QINT(K) = QI - (DELX(K) * DELQL / H) 
AINT(K) = AI ~ (DELX(K) %* DELAL / H) 
K=K +1 
GO TO 20 

END IF 


=———— CALCULATE ERROR BETWEEN ESTIMATED SLOPE AND NEW SLOPE----- 
=—<== FROM NEW INTERPOLATED VALUES----- 


E(1) = DABS(LMD(1) - (QINT(1) + AINT(1))) 
E(2) = DABS(LMD(2) - QINT(2)) 

E(3) = DABS(LMD(3) = (QINT(3) = AINT(3))) 
LMD(1) = QINT(1) + AINT(Q) 

LMD(2) = QINT(2) 

LMD(3) = QINT(3) - AINT(3) 


———— COMPARE ERROR TO ERROR TOLERANCE LEVEL, ITERATE TIL MET-----~ 
IF (CEC 1).GT.EE).OR.(E(2).GT.EE).OR.(E(3).GT.EE)) GO TO 10 


Sos CALCULATE LINEARLY INTERPOLATED VALUES OF REIMAN----- 
SSS SS VARIABLES AND MODIFIED ENTROPY AT POINT A----- 


QQINT = QQ ~ (DELX(1) * (DELQQL / H)) 
SINT = S - (DELX(2) * (DELSL 7 H)) 
RRINT = RR - (DELX(3) * (DELRRL 7 H)) 


~-=- CALCULATE SPATIAL DERIVATIVES --- 


QPRIM(1) = (QI-QIM1)/H 

QPRIM(2) = 0.DO00 

QPRIM(3) = (QI-QIM1)/H 

APRIM(1) = (AI-AIM1)/H 

APRIM(2) = 0.000 

APRIM(3) = (AI-AIM1)/H 

RETURN 

END 

SUBROUTINE COND3(QI,QIP1,AI,AIP1,RR,QQ,S,DELQQH,DELRRH,DELSH, 
Cc DELQH ,DELAH,DELT ,H,EE ,QQINT »RRINT ,SINT ,QPRIM, 
Cc APRIM,AINT ) 

382 BE RBBB BUBB BUBB BUBB BEB BUBB BEB OBBOF 

® * 

ad SUBROUTINE CONDITION 3 * 

4 x 


HEHE EI HEE IE IE IE HEHE HEHE IE IEE FETE IE IEE IE FE IE TE HEHE IE IE IE IE IE IE IE IE IE IE IE HE IE IE IE TE IE IE IE IE IE IE TEE HEE HE HEHE 


“-- FORWARD DIFFERENCE ALGORITHM FOR CALCULATING --- 
@-- RRINT,QQINT»SINT,APRIM,QPRIM,AINT ~-- 


DIMENSION LMD(3),DELX(3),QINT(3),AINT(3),E(3)},QPRIM(3 ),APRIM( 3) 
INTEGER XK 
DOUBLE PRECISION QI,QIP1,AI,AIP1,RR,QQ,S,> 

Cc DELQQH ,DELRRH ,DELSH,DELQH,DELAH, 
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Cc DELT,H,EE,AINT ,QINT,LMD,DELX,E, 
c QQINT »>RRINT >SINT >QPRIM,APRIM 


weme-TNITIAL ESTIMATE OF CHARACTERISTIC SLOPES----- 


LMO(1) = QI + AT 
LMD(2) = QI 
LMO(3) = QIP1 - AIP1 


10 K =1 

20 IF (K.LT.4) THEN 
DELX(K) = DELT * LMD(K) 
QINT(K) = QI - (DELX(K) * DELQH / H) 
AINT(K) = AI - (DOELX(K) * DELAH / H) 
K = «kK + 1 
GO TO 20 

END IF 


=---— CALCULATE ERROR BETWEEN ESTIMATED SLOPE AND NEW SLOPE----- 
a FROM NEW INTERPOLATED VALUES----- 


= DABS(LMD(1) - (QINT(1) + AINT(1))) 
E(2) = DABS(LMD(2) - QINT(2)) 


E(3) DABS(LMD(3) - (QINT(3) - AINT(3))) 
LMDO(1) = QINT(1) + AINT(1) 

LMD(2) = QINT(2) 

LMD(3) = QINT(3) - AINT(3) 


=-=<= COMPARE ERROR TO ERROR TOLERANCE LEVEL, ITERATE TIL MET----- 
IF(CE(1).GT.EE).OR.(E(2).GT.EE).OR.(E(3).GT.EE)) GO TO 10 


eccee CALCULATE LINEARLY INTERPOLATED VALUES OF REIMAN----- 
eocce VARIABLES AND MODIFIED ENTROPY AT POINT A----- 


QQINT = QQ - (DELX(1) * DELQQH 7 H) 
SINT = S - (DELKX(2) * DELSH 7 H) 
RRINT = RR - (DELX(3) * DELRRH 7/7 H) 


“-- CALCULATE SPATIAL DERIVATIVES --- 


QPRIM(1) = (QIP1-QI )/H 

QPRIM(2) = 0.000 

QPRIM(3) = (QIP1-QI )/H 

APRIM(1) = (AIP1-AI)/H 

APRIM(2) = 0.000 

APRIM(3) = (AIP1-AI)/H 

RETURN 

ENO 

SUBROUTINE CONDG(I,SHOCK ,CNTACT ,J,LNODE ,RNODE ,QQSTEP,RRSTEP,» 
C SSTEP ) 

JE IE IESE IEE JE EE HIE IEEE 36 E36 EE IEEE EIEIO F636 36 JOE 36696 IE FE FEE TOE FETE FE FE HE IEE 

% % 

* SUBROUTINE CONDITION 4 % 

% * 


HEHEHE HE IEE FE IE HE IE HEHE IE IE HE HE HE HE HEE FEI HEHE FEE HE IE IE EE HE HE IE IE HH HEE EE HE IOI ENE HE EE HE IEE 


“--- A SITUATION EXISTS AT NODE I THAT WILL BE CORRECTED IN --- 
--- SUBROUTINE CORRCT. *NODE ARRAYS ARE GIVEN INFORMATION --- 
--- ON NODE LOCATION AND SHOCK/CONTACT SURFACE PICTURE ---- 


CO4TRK - DENOTES HOW MANY NODES NEED TO BE CORRECTED 
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THIS TIME STEP 


DIMENSTON LNODE(4),RNODE(4) 
INTEGER LNODE ,RNODE ,I ,SHOCK ,CNTACT »J,CD4TRK 
DOUBLE PRECISION QQSTEP,RRSTEP ,»SSTEP 


~wo---ASSIGN FIRST NODE ENCOUNTERED DURING THE NEW TIME STEP TO ---- 
wee2=iNODE ----- 
IF(LNODE(4).LT.J) THEN 
CO4TRK = 1 
END IF 
IF(COS4TRK.EQ.1) THEN 
LNODE(1) = I 
LNODE(2) = SHOCK 
LNODE(3) = CNTACT 
LNODE(4) = J 
———— IF A SECOND NODE WITH CONDITION 4 IS ENCOUNTERED IN THE ----- 
=== SHEEP, THIS NODE IS ASSIGNED TO RNODE ----- 
EESE 
RNODE(1) = I 
RNODE(2) = SHOCK 
RNODE(3) = CNTACT 
RNODE(@) 3 J 
END IF 


aan 


QQSTEP = 0.000 
RRSTEP = 0.000 
SSTEP = 0.000 


TF(CO4TRK.EQ.1) THEN 
COGTRK = CDGTRK + 1 
ELSE 
COGTRK = 1 
END IF 
RETURN 
END 


SUBROUTINE COND5(QI ,QIM1,QIP1,AI ,AIM1L,AIP1>RR>»QQ,S,DELQRQL , 


DELQQH,DELRRL,DELRRH,DELSL ,DELSH,DELQL ,DELQH;, 
DELAL ,»DELAH,H,EE ,DELT ,QQINT ,RRINT >SINT ,AINT 5 
QPRIM,APRIM,SHOCK ,CNTACT ) 

IBBUERBUBERBRBEHBUERURB8RRBRBRBBEEBEBREBEUEBEBUERHOE 

xe * 

x SUBROUTINE CONDITION 5 x 

* % 


FE IEF IEE IEE HE HE IE FETE HEHE HEHE IEE HE IE HE IE IEIE HE IE IEIETE JE IEIE HEHE IE HEHE HEHE IE HE IEDE FE HE TE IEIE FE IE 3 IE IEE FE IE 


--- FOR SUPERSONIC FLOW WITH A DISCONTINUITY ON ONE SIDE OF THE NODE- 
~-- CALCULATES QQINT,RRINT »>SINT,APRIM,QPRIM,AINT --- 


Cc 
Cc 
C 


DIMENSION LMD(3),DELX(3),QINT(3),AINT(3),E(3),QPRIM(3),APRIM(3) 
INTEGER K,SHOCK ,;CNTACT ‘ 
DOUBLE PRECISION QI,AI,RR»QQ,S,QIM1,QIP1L,AIM1,AIP1, 
EE ,DELT,LMD,DELX,QINT ,AINT,E ,DELQL,DELQH,DELAL, 
QQINT ,SINT >RRINT »>QPRIM,APRIM,DELAH >H,» 
DELQQL ,DELQQH,DELRRL,DELRRH,DELSL ,DELSH 


--- DISCONTINUITY ON LEFT, HEADED RIGHT, NOT CROSSING NODE --- 


IF( (SHOCK .EQ.322).0R.1(CNTACT.EQ.322)) THEN 


> ITI 


aan 


ano aan aan ‘oe mem @) anoaNn 


aon 


ano 


10 K 


LMD(1) = QI ¢« AT 
LMDf2) = Qt 
LMD(3) = QIP1l - AIP1 


= 1 


20 IF (K.LT.4) THEN 


DELX(K ) 
QINT(K) 
AINT(K) 


K=K +1 
GO TO 20 


END IF 


CALCULATE ERROR BETWEEN ESTIMATED SLOPE AND NEW SLOPE 


DELT * LMDI(K) 
QI - (DELX(K) % DELQH 7 H) 
AI = (DELX(K) % DELAH / H) 


FROM NEW INTERPOLATED VALUES----- 


E(1) 
E(3) 
LMD( 1) 


LMD{ 2 ) 
LMD(3 ) 


= DABS(LMD(1) = (QINT(1) + AINT(1))) 
E(2) = DABS(LMD(2) = QINT(2)) 
= DABS(LMD(3) = (QINT(3) - AINT(3))3 


QINT(1) + AINT(1) 
QINT( 2 ) 
QINT(3) = AINT(3) 


QQINT = QQ - (DELX(1) * DELQQH / H) 
SINT = S - (DELX(2) * DELSH / H) 
RRINT = RR - (DELX(3) * DELRRH / H) 


“~~ CALCULATE SPATIAL DERIVATIVES --- 


QPRIM( 1) 
QPRIM( 2) 
QPRIM(3) 
APRIM(1) 
APRIM(2) 
APRIM( 2) 
END IF 


(QIP1-QI )/H 
0.000 
{QIP1-QI )/H 
{AIPL-AI }/H 
0.000 
(AIPI-AI )/H 


“<= DISCONTINUITY ON RIGHT, HEADED LEFT, NOT CROSSING NODE --- 


30 K 


IF( (SHOCK .EQ.232).0R.(CNTACT.EQ.232)) THEN 


INITIAL ESTIMATE OF CHARACTERISTIC SLOPES----- 


LMD(1) = QIM1 + AIM1 
LMD(2) = QI 
LMD(3) = QI - AI 


= 1 


40 IF (K.LT.4) THEN 
DELX(K) 
QINT(K) 
AINT(K) 


K 
GO 
END IF 


K + 
TO 40 
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DELT * LMDIK) 
QI - (DELX(K) * DELQL / H) 
AI - (DELX(K) %* DELAL 7/7 H) 
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awne-2-FROM NEW INTERPOLATED VALUES----- 


= DABS(LMD(1) - (QINT(1) # AINT(1))) 
E(2)} = DABS(LMD(2) = QINT(2)) 
= DABS(LMD(3} = (QINT(3) = AINT(3}))} 


= QINT(1) + AINT(1) 
LMD(2) = QINT(2) 
= QINT(3) = AINT(3) 
SS aioe COMPARE ERROR TO ERROR TOLERANCE LEVEL, ITERATE TO MEET----- 
IF ((E&01).GT.EE).OR.(E(2).GT.EE).OR.(E(3).GT.EE)) GO TO 30 


=Sa=— CALCULATE LINEARLY INTERPOLATED VALUES OF REIMAN----- 
“<< VARIABLES AND MODIFIED ENTROPY AT POINT A----- 


QQINT = QQ - (DELX(1) * (DELQQL / H)) 
SINT = S - (DELX(2) * (DELSL 7/7 H)) 
RRINT = RR - (DELA(3) * (DELRRL 7 H)) 


w-- CALCULATE SPATIAL DERIVATIVES --- 
QPRIM(1) = (QI-QIM1)/H 
QPRIM(2} = 0.D00 
QPRIM(3) = (QI-QIM1)/H 
APRIM(1) = (AI-AIM1)/H 
APRIM(2)} = 0.000 
APRIM(3) = (AI-AIM1)/H 
END IF 
RETURN 
END 
SUBROUTINE COND6( SHOCK ,CNTACT ,HALT ) 
3633 00GB BUBB EB BBB BRE EBB RR BUBB B8E 
* * 
x SUBROUTINE CONDITION 6 # 
*% % 
I3003BRE EURO BRB BEER BBR RNB BRE RE RBBB BORER REE 
INTEGER SHOCK ,CNTACT ,HALT 
PRINT * »‘THE SITUATION FOR CONDITION 6 WITH THE CONTACT SURFACE ‘ 
PRINT * ,'TO THE RIGHT OF A SHOCK BOTH HEADED RIGHT OR THE ; 
PRINT * »‘CONTACT SURFACE TO THE LEFT OF A SHOCK BOTH HEADED LEFT’ 
PRINT * ,*ADOITIONAL LOGIC IS REQUIRED TO PROCEED‘ 
PRINT * ,»‘SHOCK =*,SHOCK,* CONTACT SURFACE =',CNTACT 
HALT = 1 
RETURN 
END 
SUBROUTINE COND7( SHOCK ,CNTACT ,DELT,SIGMA ,I2,I,H»HALT>Q) 
BEBO EU BRB RBBB BUBB EE BBB RBBB BBB RBNE BBE 
* x 
% SUBROUTINE .CONDITION 7 xe 
* x 


HE IE FE FEI HE HE IE IE IE HE IE IE IE IEE HE HE JE IE IE IE IE HE IEIE HE IE IE HE FE IE HE JE IE IE HE HE HEIE IE HE IE IE IEIE IEE IE HEHE HEFCE HE IE 


DIMENSION SIGMA(4,2},I2(4) 
INTEGER HALT ,SHOCK ,CNTACT ,I2,I1 
DOUBLE PRECISION DELT,SIGMA,H>,Q 


PRINT * »'CONDITION 7 REQUIRES THAT THE CONTACT SURFACE AND : 


PRINT * ,*SHOCK(MOVING IN OPPOSITE DIRECTIONS) MEET AND CROSS. ' 
PRINT * ,»‘WHEN THEY INTERSECT,THE RESULT IS A FUNCTION OF EXIST-' 
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PRINT % »'ING CONDITICNS AROUND THEM. BOTH THE ORIGINAL SHOCK ' 
PRINT % »'AND CONTACT SURFACE WILL EXPERIENCE VELOCITY CHANGES .' 
PRINT *® »°THIS SUBROUTINE WOULD HAVE TO CALCULATE WHEN AND WHERE‘ 
PRINT %* ,»‘WITHIN THE TIME/SPACE ITERVAL THE SHOCK AND CONTACT : 
PRINT * ,‘SURFACE INTERSECT. THIS IS BASED ON KNOWN SPEED AND ' 
PRINT * ,‘THEIR RESPECTIVE LOCATIONS. THIS NEW TIME, DELT(NEN) ° 
PRINT % »'THEN COULD BE USED TO RERUN "SWEEP", WITH CONDITION 7S‘ 
PRINT * ,'EXISTING AT TIME = T + DELT(NEW). ADDITIONAL LOGIC IS’* 
PRINT % ,»*REQUIRED TO CONTINUE. ‘ 

PRINT * »'SHOCK =',SHOCK,' CONTACT SURFACE =',CNTACT 

PRINT * »‘SIGMA(1,1) =°,SIGMA(1,1),° SIGMA(1,2) =°,SIGMA(1,2) 
PRINT * »'SIGMA(2,1) =',SIGMA(2,1),° SIGMA(2,2) =',SIGMA(2,2) 
PRINT * ,'CONTACT SURFACE VELOCITY =',Q 

PRINT * ,'NODE I =°,I 

HALT = 1 

RETURN 

END 


SUBROUTINE COND/N( SHOCK ,CNTACT ,DELT »SIGMA,1I2,1,H»HALT,Q)X) 


FOIE IEIE TE EES IE IE FETE FEI IE HE HEDETE IE IE IE IE IE FEE FETE JE IEIEIE IE IE IE FEIEDE IEE IE IE JE IE IE IE IE FETE IE IE IEIEIE EE 


% 
% 
% 


% 
SUBROUTINE CONDITION 7N % 
% 


HE IED TE FE IEIE FETE IE IE IE IE IE IE IE FETE IE IE IE IE IE HE IE IE IE IE IE TE IE HE TE TE IE IE IE IE IE IETHE IE IE IE IE IE IE IE IE IE IE IE FE HE IETE 


DIMENSION SIGMA(4,2),12(4) 
INTEGER HALT ,SHOCK ,CNTACT ,I2,I 
DOUBLE PRECISION DELT,SIGMA>H,Q,xX 


PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
HALT 


RETURN 


END 


Kx KK K KK KK KK KK KK KK OK OK OX 


» ‘CONDITION 7N REQUIRES THAT WHEN EITHER A SHOCK OR ; 
» ‘CONTACT SURFACE JUMPS A NODE (AS DETERMINED BY ; 
» ‘COMPARING SIGMA(L,1) TO SIGMA(L,;2)) A CONTACT SURFACE‘ 
» OR SHOCK WOULD BE MET AND CROSSED DURING THE JUMP. : 
»' THE RESULT WHEN THEY INTERSECT IS A FUNCTION OF : 
» "EXISTING CONDITIONS AROUND THEM. BOTH THE ORIGINAL' 

» SHOCK AND THE CONTACT SURFACE WILL EXPERIENCE VELOCITY ‘ 
» CHANGES. THIS SUBROUTINE WOULD HAVE TO CALCULATE WHEN’‘ 
» WITHIN THE TIME INTERVAL AND WHERE SPACIALLY THE ? 
» "INTERSECTION OCCURS. THIS IS BASED ON KNOWN SPEED AND’ 
»' THE RESPECTIVE LOCATIONS OF THE SHOCK AND CONTACT : 
>» ‘SURFACE. THE NEW TIME, DELT(NEW}, COULD THEN BE USED ° 
»'TO RERUN "SWEEP" WITH CONDITION 4 AT THIS NODE, AND : 
» SIGHA(2,2) = SIGMA(1,2) SO THAT CONDITION 7S WOULD : 
» "RESULT IN THE NEXT TIME INTERVAL.' 

»' SHOCK =* ,SHOCK, 'CONTACT SURFACE =' ,CNTACT 

» ‘SIGMA(1,1) =*°,SIGMA(1,1),° SIGMA(1,2) =',SIGMA(1,2) 

» SIGMA(2,1) =',SIGMA(2,1),°' SIGMA(2,2) =* ,SIGMA(2,2) 

» CONTACT SURFACE VELOCITY =',Q 

»' NODE I =',I,* LOCATION AT X =' 4X 


SUBROUTINE COND7S( SIGMA ,HALT ,SHOCK ,CNTACT ) 


FEIEH HE HEIEHE HE HEHE IE HE IE IE HEHEHE HEHEHE HEHE IE EIE HE HE HE ETE HE HE HE HEE HE HEHEHE HEHEHE IE HEHE HEHE EEE HE HEHE HHH 


% 
% 
% 


% 


SUBROUTINE CONDITION 7S % 
x 


HEHE IE IE IE HE HE IE HE HE IE IE IE HE IEE HE IE IE HE HE IE IE IE IE IE IE IE IE IE IEE HE HE HE IE HE HE IETE HE HEE FETE IE FEI HEHE HEHE IE HEHEHE 


DIMENSION SIGMA( 4,2) 
INTEGER SHOCK ,CNTACT ,HALT 


DOUBLE 


PRINT % ,*CONDITION 7S HAS BEEN MET. 


PRECISION SIGMA 


THIS MEANS THAT THE SHOCK '* 


PRINT % »"AND CONTACT SURFACE ARE LOCATED AT THE SAME X AT A : 
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PRINT * »‘'TIME OTHER THAN ZERO. THIS SUBRGUTINE WOULD HAVE TO ' 
PRINT ® ,'DETERMINE THE RESULT OF THE INTERSECTION BASED ON THE ° 
PRINT ® »‘CONDITIONS TO THE LEFT AND RIGHT. ADDITIONAL LOGIC IS' 
PRINT # ,*REQUIRED TO PROCEED. ' 

PRINT * »‘SHOCK =',SHOCK,* CONTACT SURFACE =‘ ,CNTACT 

PRINT * »‘SIGMA(1,1) =',SIGMA(1,1),° SiGMA(2,1) =',SIGMA(2>1) 
HALT = 1 

RETURN 

END 


SUBROUTINE COND8{ SHOCK ,CNTACT ,HALT } 


HEFCE HE FEE IE HE HEI IE IE IE IE IE IE IE IE IE HEI IE IEE HEHE HEE IE IE IE IE IE IE IE IE IH EI HE IH IH HIE HHH IE HE HHH HEHE 


* * 
* SUBROUTINE CONDITION 8 * 
& % 


FE IEIE FE IE HE ETE IE HE IE HIE HE HE IEE HE HEHE EIEE TE EF IE IE IE HE HE FEE IE HE IE HEI IE HE HE IE TEE HE IE IEE EE HEE HEE HH 


INTEGER HALT »SHOCK ,CNTACT 


PRINT * ,»,'CONDITION 8 COULD RESULT ONLY AFTER THE ORIGINAL SHOCK‘ 

PRINT * »'HAS CROSSED THE CONTACT SURFACE, AND THE SUBSEQUENT : 

PRINT * »'CONDITIONS DETERMINED. ADDITIONAL LOGIC IS REQUIRED TO' 

PRINT * » ‘CONTINUE. ' 

PRINT * ,»‘'SHOCK =',SHOCK,' CONTACT SURFACE =',CNTACT 

HALT = 1 

RE TURN 

END 

SUBROUTINE CORRCT(LNODE ,RNODE »N,SIGMA,H,QQ,RR»S5G,G1,G62,12,X2.W> 
Cc AR »,DQ,VS,A,Q) 

BE HE EE SEE EE he RG EI RE ES GE SE IE EE EI I A SEED 

% *% 

* DISCONTINUITY CORRECTION SUBROUTINE % 

~ % 


JEIEIE HE HE IE IE HE IE IE IE HE HE HE IE IE HE IE IE FEE IE HE HE IESE IEE IE EE HEE HE IE IE EE IE IE HE EE HEE HE HEHEHE HHH HH 


~----------- VARIABLE DEFINITIONS ------------ 


ENTRPY = MODIFIED ENTROPY 

SNDSPD - SONIC VELOCITY 

UA = VELOCITY RELATIVE TO THE SHOCK, LEFT SIDE 
UB = VELOCITY RELATIVE TO THE SHOCK, RIGHT SIDE 
VLCTY = VELOCITY 


INTEGER LNODE ,RNODE »N»I2,NODE »SHOCK »CNTACT 5K 
DIMENSION LNODE(4),RNODE(4),SIGMA(4,2),QQ(N),RRIN)»,S(N)J,I2(4), 
Cc X2(4),A(N),QOUN) »XA(4),XB(G) 
DOUBLE PRECISION SIGMA,QQ,RR>»S>A,Q, 
H,G,G1,62,X2,N,AR,0Q,VS, 
RRA,RRB »,QGA »QQB,AA,AB»,QA,QB,SA»SB> 
SA1,SAZ,VLCTY »SNOSPD ,ENTRPY »QQCALC »RRCALC > 
XB XA 


Na O 


QQCALC(VLCTY ,SNDSPD,ENTRPY) 2 VLCTY + SNDSPD*ENTRPY 
RRCALC(VLCTY »SNDSPD ,ENTRPY) = VLCTY =- SNDSPDXENTRPY 


RRA = 0.000 
RRB = 0.000 

QQ@A = 0.000 

Q@B = 0.000 

AA = 0.000 

; AB = 0.000 
QA = 0.000 


as 


aaAaAno 


aan 


5 


0.000 
0.000 
0.000 
DO 15 K=15% 


w 
> 
to 0 Ot 


AA(K}=0 .00 
XB(K 3=0.00 


CONTINUE 


“NODES. ALSO, IF SHOCK AND CONTACT SURFACES ARE CLOSE----- 
“ENOUGH(WITH IN 2H) TO INTERACT----- 


IF( ( LNODE(2).EQ.100).0OR.{LNODE(3).EQ.100)) THEN 
IF(( LNODE(1).GE.(RNODE(1)-2)).AND.(RNODE(1).GT.0)} THEN 


ELSe 


ENO IF 


GO TO 20 


GO TO 30 


ELSE IF(RNODE(1).GT.0O) THEN 
GO TO 20 


ELSE 


GO TO 10 


END IF 


10 SHOCK = LNODE( 2) 
CNTACT = LNODE(3) 


---SH 


OCK ON LEFT, HEADED RIGHT, JUMPS OR NOTSCONTACT SURFACE ON----- 
---RIGHT, HEADED LEFT, DOESN'T CROSS NODE --- 


IF( (SHOCK.EQ.322.0R.SHOCK .EQ.321).AND.(CNTACT.EQ.232)} THEN 


I2(1) 
I2(2) 


LNODE(1) + 1 
LNODE(1) + 1 


CALL DELTAX(1I2,SIGMA,XB,XA,H) 


IF (12(2).EQ.N) THEN 


ELSE 


CALL 8BORY(RR(N),QQ(N),S(N),RRB,QQB,SB,AB,QB) 


CALL EXTRAP(RR(12(023),RRIT20 2 34+13,QQ(T2(2)),QQ(IT2(2)4+1), 


END IF 


S(1T20233,S012(2)3+1),xXB(2),H,RRB»QQB,SB,AB,QB) 


SA = S(I2( 23-1) 
CALL CSJUMP(AB,QB,SB,SA,G2,QA,AA ) 


IF SHOCK INTERACTS CALCULATE CHANGE IN VARIABLES ACROSS SHOCK---- 


IF(SHOCK.EQ.321) THEN 


QB = QA 

AB AA 

$B SA 

CALL SKJUMP(AB,Q@B,SB,AR,0Q,VS,G,G1,W,AA,QA,SA) 
QQA = QQCALC(QA,AA,SA) 

RRA = RRCALC(QA,AA,SA) 


IF (12(1)}.EQ.2) THEN 
CALL BBDRY(RRA,QQA,SA,RRIT2(1)-1),QQ(T2(1)-1)3, 
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Cc S(T201 )-1),A0T201)-1),Q(0T2(1)-1)) 

ELSE 

CALL INTERP( RRA,RR(T2(1 )-2),QQA,QQ(T2(1)-2),SA, 
S(T201)-2),XAC1),0XAC1 4H),RRIT201)-21), 
QQ(T2(1 )-1),S0T201)-1),A0T201)-1), 
Q(IT2(1)-1)) 

END IF 


aan 


ELSE 
QQ(I2(2)-1) 
RR(I2( 2 )-1) 
S(I2(2)-1) 
A(T2(2)-1) 
Q(IT202)-1) 


QQCALC(QA,AA,SA) 
RRCALC(QA,AA »SA) 
SA 
AA 
QA 


END IF 


---SHOCK ON RIGHT, HEADED LEFT, CROSSES NODE OR NOT; CONTACT SURFACE- 
~--ON LEFT,HEADED RIGHT, DOES'NT CROSS NODE--- 


ELSE IF((SHOCK.EQ.232.OR.SHOCK.EQ.231).AND.(CNTACT.EQ.322)) THEN 


aan ana 


T2(1) LNODE(1) 
T2(2) LNODE(1) 
CALL DELTAX(I2 ,SIGMA »XB >XA 5H ) 


IF (I2(2).EQ.2) THEN 
CALL BBDRY(RR(1),QQ(1)5,501)>RRA,QQA,SA,AA,QA ) 
ELSE 
CALL EXTRAP(RR(I2(2)-1),RRIT21 2 J-2) ,QQ(T2( 2)-1),QQ(T2(2)-2), 
Cc SCT202 J-1),S50T 202 )-2),XA02)5H»RRAQQA,SA,AA,QA ) 
END IF 


SA = S(T202)) 
CALL CSJUMP( AA ,QA,SA>5SB ,G2,QB,AB ) 


—————— IF SHOCK INTERACTS CALCULATE CHANGE IN VARIABLES ACROSS SHOCK---- 


IF(SHOCK .EQ.231) THEN 
QA = GB 
AA = AB 
= SB 
CALL SKJUMP( AA,QA,SA,AR,DQ,VS»G,G1>N,AB ,QB »SB ) 
QQB = QQCALC(QB,AB,S8 ) 
RRB = RRCALC{ QB,AB,SB3 ) 


IF (I2(1).EQ.N) THEN 
CALL BBORY( RRB,QQB,SB,RR(T2(1)),QQ(T2(1)),S5(T2(1)); 
Cc A(T2(1)),Q(I2¢1))) 
ELSE 
CALL INTERP(RRB,RRII211)41),QQB,QQ(I2(1)4+1),SB8, 
C SCT2C1 )+1),XB(1),0X801)4+H) »RRIT201)),; 
Cc QQ(T2(1)),S(T201)),A0T201)),Q0T201))) 
END IF 
ELSE 
QQ(T2(2)) = QQCALC(QB,AB,SB) 
RR(I2(2)) = RRCALC(QB,AB,SB) 
S(T202)) SB 
A(T2(2)) AB 
Q(T212)) QB 


END IF 
C 


C --- SHOCK ON RIGHT,HEADED RIGHT, DOESN'T CROSS OR ON LEFT ,HEADED--- 
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C --- RIGHT, DOES CROSS; AND CONTACT ON LEFT,HEADED RIGHT, CROSSES--- 


C --- GR NOT: OR SHOCK ON RIGHT, HEADED LEFT, DOESN'T CROSS NCCE --- 
C --- WITH CONTACT ON LEFT, HEADED RIGHT, AND CROSSED NCDE --- 
C 
ELSE IF(( (SHOCK.EQ.222.0OR.SHOCK.EQ.321).AND.(CNTACT.EQ.321.0OR. 
Cc CNTACT.EQ.322 )).0R.(SHOCK.EQ.232.AND.CNTACT.EQ.321)) 
Cc THEN 
Cc 
C -r--- DETERMINE NODE TO RIGHT OF SHOCK AND CONTACT SURFACE----- 
C 
T2(¢1) = LNODE(1) + 1 
T2(2) = LNODE(1) + 1 
CALL DELTAX(I2,SIGMA ,XB,xXA>H) 
Cc 
C ----- EXTRAPOLATE TO RIGHT FACE OF SHOCK----- 
Cc 
IF (3201).EQ.N) THEN : 
CALL BBDRY(RRIN),QQ(N),S(N),RRB,QQB,SB,AB ,>&B ) 
ELSE 
CALL EXTRAPCRR(T2(1)),RROIT201)41),QQ(T2(1)),QQ(T2(1 +1), 
Cc $(T201)),S0T201)+41),XB01),H,RRB,QQB,SB8 ,AB ,QB ) 
END IF 
Cc 
C e---- CALCULATE CHANGE IN VARIABLES ACROSS SHOCK----- 
Cc 
IF(SHOCK.EQ.232) THEN 
AA = AB/AR 
SAl1l = (G1/G)*DLOG( (2.D00*G( W*¥2 )-G+1.D00) 7 (G+1.D000)) 
SA2 = GI*¥DLOG(((G-1.D00 )¥(W*¥*2)42.000) 7 ((G+1.D00 )x 
Cc (W*#2))) 
SA = SB4SA1+SA2 
UB = QB - VS 
UA = UB -AA*DQ 
QA = UA + VS 
ELSE 
CALL SKJUMP(AB,QB,SB>,AR,DQ,VS>G,G1>N,AA ,QA,SA ) 
END IF 
C 
C ee-e- CALCULATE CHANGE IN VARIABLES ACROSS CONTACT SURFACE----- 
C 
IF(CNTACT.EQ.322) THEN 
QQITZ(1lj-1) = QQCALC(QA,AA,SA) 
RRII2( 1-1) = RRCALC(QA,AA,SA) 
S(T2(1l)J-1) = SA 
A(I2(1)-1) = AA 
QiT2(1ljJ-1) = QA 
ELSE 
QB = QA 
AB = AA 
$B = SA 
SA = S(I2(2)-2) 
CALL CSJUMP( AB ,QB,SB,SA,G2,QA,AA ) 
QQA = QQCALCIQA,AA,SA) 
RRA = RRCALC(QA,AA,SA) 
¢ 
C o---- INTERPOLATE TO NODE(IT2(2)-1) AND ASSIGN CORRECTED VALUES----- 
Cc 
IF (12(2).EQ.2) THEN 
CALL BBDRY( RRA ,QQA,SA,RRIT2( 2 J-1) ,QQ(I2(2)-1), 
Cc S$(T202 )-1) ,ACT202 )-1),Q(0I202)-1)) 
ELSE 
CALL INTERP(RRA,RR(I2( 2 )-2) ,QQA ,QQIT2(2)-2),S5A, 
Cc S$(I202 )-2),XA(2),0XAC2 4H) ,RRIT202 )-1)> 
Cc QQ(T2(2)-1),S0T202)-1),A0T202)-1), 
Cc Q(I202)-1)) 
END IF 
END IF 
Cc 


C ---SHOCK ON LEFT, HEADED LEFT, DOESN'T CROSS NODE, OR ON RIGHT >--- 
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C ---HEADED LEFT, CROSSED NODE$} AND CONTACT GN RIGHT, HEADED LEFT,---. 
C ---CROSSES NODE OR NOT: OR SHOCK ON LEFT, HEADED RIGHT, DCESN'T--- 
C ---CROSS NODE, AND CONTACT ON RIGHT, HEADED LEFT, CROSSED NODE--- 
Cc 
ELSE IF(( (SHOCK .EQ.332.0R.SHOCK.EQ.231). AND. (CNTACT.EQ.232.0R. 


C CNTACT.EQ.231)).0OR. (SHOCK .EQ.322.AND.CNTACT.EQ.231)) 
C THEN 
Cc 
C wr--- DETERMINE NODE TO RIGHT OF SHOCK AND CONTACT SURFACE----- 
C 
T2(1) = LNODE(1) 
I2(2) = LNODE(1) 
CALL DELTAX(I2,SIGMA,XB>,XA>H ) 
S 
C e----EXTRAPOLATE TO LEFT FACE OF SHOCK----- 
C 
IF (I12(1).EQ.2) THEN 
CALL BBDRY(RR(1),QQ(1),S(01),RRA»,QQA,SA,AA ,QA) 
ELSE 
CALL EXTRAP(RR(I2(1 )-1),RR(I2(1)-2),QQ(I2(1)-1) ,QQ(IT2(1)-2), 
C $(I201 )-1),S0I201 J-2),XA0C1)5H,RRA,QQA,SA,AA,QA) 
END IF 
C 
C -----CALCULATE CHANGE IN VARIABLES ACROSS SHOCK----- 
Cc 
IF (SHOCK .EQ.322) THEN 
AB = AA/AR 
SA1 = (61/G )¥DLOG( (2. D00xG( W*x*2 )-G+1.000) / (G+1.D0D00)) 
SA2 = G1¥DLOG( ((G-1.D00 )*¥( Wx*2)4+2.000) 7 (1(64+1.000)% 
od { W332) )) 
SB = SA+SA1+4SA2 
UA = QA - VS 
UB = UA - AB*DQ 
QB = UB + VS 
ELSE 
CALL SKJUMP(AA,QA,SA,AR,DQ,VS,G,G1,W,AB >QB,SB ) 
END IF 
Cc 
C +---- CALCULATE CHANGE IN VARIABLES ACROSS CONTACT SURFACE----- 
C 
IF(CNTACT.EQ.232) THEN 
QQ(I2(1)) = QQCALC{ QB,AB,SB) 
RR(I2(1))} = RRCALC(QB,AB,SB) 
$(I2¢1)) = SB 
A(I2(1)) = AB 
Q(I2(1)) = QB 
ELSE 
QA = QB 
‘AA = AB 
SA = SB 
SB = S$(I2(2)41) 
CALL CSJUMP( AA,QA,SA,SB,G62,QB,AB ) 
QQ@B = QQCALC(QB,AB,SB ) 
RRB = RRCALC(QB,AB,SB) 
c 
C ----- INTERPOLATE TO NODE(I2(2)) AND ASSIGN CORRECTED VALUES----- 
c 
IF (I12(2).EQ.N) THEN 
, CALL BBDORY( RRB ,QQB ,SB ,RR(I2(2)),QQ(I2(2)),S(I202)); 
C A(I2(2)),Q(I202))) 
ELSE 
CALL INTERP(RRB,RR(I2{ 23}4+1),QQB,QQ(I2(2)+1),SB, 
C S(I2(02 )41),xXBC2 ),0XB(2 14H), ,RRII2(2)), 
Cc QQ(I2(2)),S(1202)),A0T202 3),Q0IT202))) 
END IF 
END IF 
END IF 
GO TO 40 . 
Cc 
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LEFT NODE: 

m--SHOCK ON RIGHT, HEADED RIGHT, JUMPS NODE WITH CONTACT ON LEFT--- 
---HEADED RIGHT, JUMPS NODE CR NOTs OR NO SHOCK WITH CONTACT ON --- 
---LEFT, HEADED RIGHT, JUMPS NODE--- 

RIGHT NODE: 

“--SHOCK ON LEFT, HEADED RIGHT, JUMPS NODE WITH NO CONTACT SURFACE--- 


20 IFC (( LNODE(2).EQ.221).AND.(LNODE(3).EQ.321.0R.LNODE(3).EQ.322)) 
Cc .OR.( LNODE(2).EQ.100.AND.LNODE(3).EQ.321)) THEN 
IF ( (RNODE(2).EQ.321).ANOD.(RNODE(3).EQ.100)) THEN 


I2(1) = RNODE(1) + 1 
IF( LNODE(2)}.EQ.322) THEN 
I2(2) = LNODE(1) 
ELSE 
I2(2) = LNODE(1) +1 
END IF 
CALL DELTAX(I2,SIGMA »>XB,XA,H) 


IF (32(1).EQ.N) THEN 
CALL BBORY(RRIN),QQ(N),S(N),RRB»,QQB,SB,AB,GB) 


ELSE 

CALL EXTRAP(RR(IT201)),RR0I1T201)4+1),QQ(T2(1)),QQ(T2(1)+1);, 
C S(IT2013),S0I201)+1),4B(1),H,RRB,QQB,SB,AB ,QB ) 

END IF 


IF ( LNODE(2).EQ.100) THEN 
QQA = QQCALC(QA,AA,SA) 
RRA = RRCALC(QA,AA,SA) 


=o2>= INTERPOLATE TO NODE(I2(1)-1) AND ASSIGN CORRECTED VALUES----- 


IF (12(1).EQ.2) THEN 
CALL BBDRY(RRA,QQA,SA,RRIT2(01 J-1),QQ(I2(1)-1), 
C S(T201)-1),A0T201 )-1),Q0¢T201)-1)) 
: ELSE 
CALL INTERP(RRA,RRIT2(1)-2),QQA,QQ(T2(1)-2),SA, 
S¢IT201)-2),XA(1),0XA01)+H)I,RRIT201)-1), 
QQ(T2(1)-1),S0¢T201)-1),A0T2(01)-1), 
Q(T201)-1)) 
END IF 


aA 


IF (T2(2).EQ.N) THEN 
CALL BBDRY(RRIN),QQ(N),S(N),RRB»QQB,SB,AB ,QB ) 


ELSE 
CALL EXTRAP(RR(T2(2)),RRIT202)4+1),QQ(T2(2)),QQ(T2(2)+4+1), 
C S(T202)),S01T202 )+1),XB(2),H,RRB,QQB,SB,AB,GB ) 
END IF 
ELSE 
QQ(T2(1)-1) = QQCALCIQA,AA,SA} 
RR(I2(1L)-1) = RRCALC(IQA,AA,SA} 
S(I2(1l)J-1) = SA 
A(T2(1l)J-1) = AA , 
Q(T2(1)-1) = QA 
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IF( LNOOE(3).EQ.322) THEN 
QQ(T2(2)) = QQrl2(1)-1) 
RR(IT2(2)) = RR(T201)-1) 

S(Tet2)) = StT201)-1) 
A(I2(2)) A(I2(1)-1) 
Q(T202)) Q(IT2(1)-1) 
GO TO 40 


ELSE 
QB 
AB 
$B 


a uf 8 
> 
> 


ENO IF 
END IF 
SA = S(I2(2)-2) 
CALL CSJUMP( AB ,QB,SB,SA,G2,QA>,AA ) 
QQA = QQCALC(IQA,AA,SA) 
RRA = RRCALC(QA,AA,SA) 
Cc 
C -----INTERPOLATE TO NODE(I2(2)-1) AND ASSIGN CORRECTED VALUES----- 
S 
IF (31202).EQ.2) THEN 
CALL BBORY(RRA,QQA,SA,RRII2( 2 J-1) ,QQUI2(2 J-1); 
= S(I202 )-1),A0T202 J-1),Q0I2(2 J-1)) 
ELSE 
CALL INTERP(RRA,RRII2( 2 )-2 ) ,QQA,QQ(I212)-2),SA, 
S$(I202 3-2) ,XA02),0XAC2 4H) RRIT202 0-1), 
QQ I2(2)-1),S0I202 J-1) ,ACI202)-1), 
Q(I2(2)-1)) 
ENO IF 


aan 


ENO IF 


LEFT NODE: 
---SHOCK ON LEFT, HEADED RIGHT, JUMPS NODE WITH NO CONTACT--- 
RIGHT NODE: 
“--NO SHOCK, WITH CONTACT ON RIGHT, HEADED LEFT ,JUMPS NODE~--- 


AAaAAAO 


ELSE IF((LNODE(2).EQ.321).AND.(LNODE(3).EQ.100)) THEN 
IF(( RNODE(23.EQ.100).ANOD.(RNODE(3).EQ.231)) THEN 


wm~-=DETERMINE NODE TO RIGHT OF SHOCK AND CONTACT SURFACE----- 


aan 


EZCY) LNODE(1) + 1 
T2(2) RNODE(1) 
CALL DELTAX(1I2,SIGMA »XB>XA>H ) 


IF(I201).EQ.(1202)-1)) THEN 
AA A(T201)) 
Q(T201)) 
S(T201)) 
S(T202)4+1) 
CALL CSJUMP(AA,QA,SA,SB,G2,QB, AB ) 
QQB = QQCALC(QB,AB,SB) 
RRB = RRCALC(QB,AB,SB) 


YW 
> 
la ott ot 


IF (12(2).EQ.N) THEN 
CALL BBORY( RRB ,QGB,SB,RRIT2(2)),QQ(T202)),S(1202)); 

Cc ACT202)),Q(I202))) 

ELSE 

CALL INTERP(RRB,RRIT2( 2 +1), QQB,QQ(I2( 2 )+1),S8, 
Cc SOI202)41),XB(2 ),(XB(2 4H)» RRIT202)); 
Cc QQ(12(02)),S01202)),A0T202)),Q0T202))) 

END IF 
A(T2(1)) 
Q(I2(01)) 
S¢I201)) 


&) 
w 
| ie | ee 
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CALL SKJUMP( AB ,QB,SB,AR,0Q,VS,G,G1,W,AA,QA,SA) 
QQA ® QQCALC(GQA,AA>SA) 
RRA = RRCALC(QA,AA,SA) 


IF (12(1).EQ.2) THEN 
CALL BBORY( RRA, QQA,SA,RRIT2(1)-1),QQ(T2(1)-1), 
C S(I2¢1 )-1) ,ACT201 J-1),Q(T201 )-1)) 
ELSE 
CALL INTERP(RRA,RR(1T2(1)-2) ,QQA,QQ(T2(1)-2),SA, 
S(¢T201)-2),XA(13,0XA0C1I4H),RRIT201)-1), 
QQ(T2(1)-1),S¢T201l)-1),A¢T2¢1)-1), 
Q(T2(1)-1)) 
ENO IF 


aan 


ELSE 

RRIT2( 1-1) = RRIT2(1)-2) 

QQ(T201)-1) = QQi(T2(1)-2) 
S(T201)-1) S(I2(1)-2) 
A(T2(1)-1) A(T2(1)-2) 
Q(T2(1)-1) Q(I2(1)-2) 

RRII2(2)) = RRIT202)41) 

QQ(T202)) = QQiT2(2)+1) 
 ~ObZt2)) si tzt2)+l) 
A(T2(2)) A(I2(2)+1) 
Q(I2(2)) QiT2(2)+1) 


ENO IF 
END IF 


RIGHT NODE: 

---SHOCK ON RIGHT, HEADED LEFT, JUMPS NODE WITH NO CONTACT--- 
LEFT NOOE: 

“--NO SHOCK, WITH CONTACT ON LEFT, HEADED RIGHT, JUMPS NODE--- 


ELSE IF((RNODE(2).EQ@.231).AN0.(RNODE(33.EQ@.1003) THEN 
IF(( LNODE(2).EQ.100).AND.{LNODE(3).EQ.321)) THEN 


I2(1) = RNODE(1) 
I2(2) = tNODE(1) + 1 
CALL DELTAX(I2,SIGMA,xXB,>XA>H) 


AB = A(T2(2)) 

QB = Q(T2(2)) 

SB = S(T202)) 

SA = S(T2{2)-2) 

CALL CSJUMP( AB ,QB,SB,SA,G2,QA,AA) 
QQA = QQCALC(QA,AA,SA) 

RRA = RRCALC(QA,AA,SA) 


IF (I12(2).EQ.2) THEN 
CALL BBDRY! RRA, QQA,SA,RRIT2( 2)-1),QQ(T2(2)-1), 

C S(T202)-1),A0T2(2)-1),Q(T2(2)-1)) 
ELSE 
CALL INTERP(RRA,RR(IT2( 2 3-2) ,QQA,QQ(T 21 2)-2),SA, 

S¢I202 322 ),XA023,0XAC2 4H) RRII202 3-1), 

QQIT202)-1),S01202 J-1) ACT 202-1); 

Q(I2(2)-1)) 


Noa 0 


ENO IF 

A(T2(2)) 

Q(T2(2)) 

S(I2(02)) 

CALL SKJUMP(AA,QA,SA,AR,0Q,VS,G,G1,W,AB,QB,SB ) 


) 
> 
Te 
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ocoee INTERPOLATE TO NODE(I2(1)) AND ASSIGN CORRECTED VALUES----=- 
IF (I2(01).EQ.N} THEN 
CALL BBDRY(RRB,QQB,SB,RRIT2(1}),QQ(I2(1)),S(T2(1)), 
Cc ACI2(01)),Q(T2013)) 
ELSE 
CALL INTERP(RRB,RR(I2(13+1),QQB,QQ(12(1)4+1),S8, 
C S$(I2(1)+13,xXB(1),(0XB01)4H),RRII201)3, 
C QQ(I2(01)3,S(T2013),A0T201)3,Q0T201)3) 
END IF 
ELSE 
RR(IT2(2)-1)} = RRIT2(2 3-2) 
QQ(I2(2}-13 = QQ(T2( 2 3-2} 
S(I2(2}-1}3 = S(T2(2)-23 
A(T2(2)-1) = A(I2(2}-2) 
Q(T20 23-1) = Q(T2( 23-2) 
RR(I2(13}}3 = RRIT2(1 341) 
QQ(T2(1)3} = QQ(T2(1)+1} 
StI2(1})} = S(T201)4+13 
A(I2(1)}3 = A(I2(1}3413 
Q(T2(1)) = QtIT2(13+1) 
END IF 
END IF 
RIGHT NODE: 
---SHOCK ON LEFT, HEADED LEFT, JUMPS NODE WITH CONTACT ON RIGHT-=-- 
@---HEADED LEFT, DOES NOT CROSS NODE OR NO SHOCK WITH CONTACT ON--- 
---RIGHT, HEADED LEFT, CROSSED NODE--- 
LEFT NODE: 
---SHOCK ON RIGHT, HEADED LEFT, JUMPS NODE WITH NO CONTACT--- 
ELSE IF(({RNODE(2).EQ@.3313.AND.0RNOOE(3).EQ.231.0R.RNODE(33}.EQ. 
C 232)).OR.(RNODE( 2).EQ.100.AND.RNODE(33.EQ@.231)3 THEN 
IF({ LNODE(2).EQ.231}3.AND.{LNODE(3}.EQ.100}3} THEN 
2ecee DETERMINE NODE TO RIGHT OF SHOCK AND CONTACT SURFACE----- 
T2(1) = LNODE(1}3 
IF( RNODE(3).EQ.232)} THEN 
I2(2} = RNODE(1)} + 1 
ELSE 
I2(2} = RNODE(1} 
END IF 
CALL DELTAX(I2,SIGMA,xXB,XA,H} 
ocece EXTRAPOLATE TO LEFT FACE OF SHOCK----- 
IF (12(13.EQ.2) THEN 
CALL BBDRY(RRI1),QQ(1),S01},RRA,QQA,SA,AA,QA ) 
ELSE 
CALL EXTRAP(RR(IIT2(1)J-1),RROT2(1)-2) QQ(I2(1i-1},QQ(I2( 13-2), 
Cc S(I2(13-1),S0I201 J-23,XAC1),H»,RRA,QQA,SA,AA>GQA } 
END IF 
ceere CALCULATE VARIABLE CHANGE ACROSS SHOCK----- 
CALL SKJUMP(AA,QA,SA,AR,DQ,VS,G,G1,W,AB,QB,SB) 
-o2-- DETERMINE CORRECTED VALUES AT NODE(I2(1}}3 AND ASSIGN THEM----- 
IF(RNODE(2).EQ.100} THEN 
QQ@B = QQCALC{(QB,AB,SB} 
RRB = RRCALC(QB,AB,SB) 
e2ere INTERPOLATE TO NODE(I2(1)) AND ASSIGN CORRECTED VALUES----- 


QGB = QQCALC(6B,A8B,SB) 
RRB = RRCALC(QB8,A3,SB ) 
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IF (242(1).EQ.N) THEN 


CALL BBDRY(RRB,GQB,SB,RR(IT2(1)3,QQ(I2(1)),S(I2(1)), 
c A(T2(1)),Q(T2(1))) 
ELSE 
CALL INTERP(RRB,RR(T2(1)41),QGB ,QQ( 72(1)41),SB, 
Cc $(I201)41),XB(1),0XBC1)4H) »RRII201)), 
Cc QQ(T2(1)),S(T201)),A0T201)),Q0T2(01))) 
END IF 
[ 
C e---- EXTRAPOLATE TO LEFT FACE OF CONTACT SURFACE----- 
C 
IF (132(2).EQ.2) THEN 
CALL BBDRY(RR(1),QQ(1),S(1),RRA,QQA,SA,AA,QA) 
ELSE 
CALL EXTRAP(RR(I2( 2 )-1),RR(T2( 2 J-2 ) ,QQ(T2(2)-1),QQ(I2(2)-2), 
C $(T2(2)-1),S¢I202 J-2),KA( 2 3,H,RRA»,QQA,SA,AA QA) 
END IF 
ELSE 
QQ(I2(1)) = QQCALC(QB,AB,SB) 
RR(I2(1)) = RRCALC(QB,AB,SB) 
$(I2(1)) = SB 
A(I2(1)) = AB 
Q(T2(1)) = QB 
IF(RNODE(3).EQ.232) THEN 
RR(IZ2(1)41) = RRII201)) 
QQ(I2(1)4+1) = QQ(I2(1)) 
$(I2(1)41) = S$(I2(1)) 
A(T2(1)41) = A(I2(1)) 
Q(I2(1)41) = Q(I2(1)) 
GO TO 40 
ELSE 
QA = QB 
AA = AB 
SA = SB 
END IF 
ENO IF 
$B = $(IT2(2)+41) 
c 
C ser-- CALCULATE VARIABLE CHANGE ACROSS CONTACT SURFACE----- 
C 
CALL CSJUMP(AA,QA,SA,S8,G2,QB,AB ) 
QQB = QQCALC( QB,AB,SB} 
RRB = RRCALC(QB,AB,SB) 
c 
C seer" INTERPOLATE TO NODE(I2(2)) AND ASSIGN CORRECTED VALUES----- 
C 
IF (I12(2).EQ.N) THEN 
CALL BBDORY(RRB,QQB,SB,RR(I2(2)),QQ(IT2(2)),S(I2(2)), 
Cc AC(I2(2)3,Q(T2(2))) 
ELSE 
CALL INTERP(RRB,RRII2( 2 )41),Q8Q8,QQ(I2(2)4+1),S8, 
C $(I2(2)41),XB(2),(XB(2 )4H) »,RRIT202)); 
C QQ(I2(2)),S(1T2(2)),A0T2023),Q(T202))) 
END IF 
END IF 
END IF 
GO TO 40 
c 
C --- BRANCH HERE IF THERE ARE ONE OR TNO NODES TO CORRECT, BUT THEY --- 
C --=- ARE SEPERATED BY MORE THAN 2H. CHECK FOR SHOCK/CONTACT SURFACE--- 
C --- INTERACTION, AND CORRECT APPROPRIATELY --- 
Cc 


aan 


30 NODE = LNODE(1) 
SHOCK = LNODE( 2) 
CNTACT = LNODE(3) 


--- NO SHOCK, CONTACT SURFACE ON LEFT, HEADED RIGHT, JUMPS NODE --- 


35 IF((SHOCK.EQ.100).AND.(CNTACT.EQ.321)) THEN 
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aa 


I2(2) = NODE + i 
CALL DELTAX(I2,SIGIHA,XB,XA>H ) 


wem--CHECK FOR A SHOCK INTERACTING WITH CONTACT SURFACE AT NODE----- 


IFC (SIGMAC1,2).LT.(X202 34H) ). AND. 


Cc (SIGMA(1,2).S6T.(X2(2))) JTHEN 
AB = A(T2(2)) 
QB = Q(I2(2)) 
SB = S$(I2(2)) 
ELSE 


IF (12(2).EQ.N) THEN 
CALL BBORY(RR(N),QQ(N),S(N),RRB,QQB,SB,AB,GB ) 
ELSE 
CALL EXTRAP(RR(T202)),RRIT202 341), QQ(T202)),QQ(T2(2)41), 
Cc $(I202)),S(1T202 )+1),X802),H,»RRB,QQB,SB,AB,GB ) 
END IF 
END IF 
SA = S$(I2(2)-2) 


\/-- CALCULATE VARIABLE CHANGE OVER CONTACT SURFACE --- 
CALL CSJUMP(AB,QB»SB,SA,G2»QA,AA ) 


QQA = QQCALC(QA,AA,SA) 
RRA RRCALC(QA,AA,SA) 


IF (1202).EQ.2) THEN 
CALL BBORY(RRA,QQA,SA,RRII202 )-1) ,QQ(T2(2)-1), 
c S(I202)-1),A0T202 )-1),Q(T202)-1)) 
ELSE 
CALL INTERP(RRA,RRIT2( 2 J-2),QQA,QQ(T2( 2 J-2),5A, 
S¢T202)-2) ,XA02),0XA0 2 4H) ,RROIZ202 3-1), 
QQ(T2(2)-1),50T202 J-1),A0T202)-1), 
Q(IT2(2)-1)) 


aaa 


ENO IF 
--- NO SHOCK, CONTACT SURFACE ON LEFT, HEADED LEFT, JUMPS NODE --- 
ELSE IF( (SHOCK.EQ.100).AND.(CNTACT.EQ.231)) THEN 


I2(2) = NODE 
CALL DELTAX(I2,SIGMA,XB,XA,H) 


IFC CSIGMA(1,2).GT.0X202 J-(2.000%H) }).AND.(SIGMAC1,2).LT. 


c (X202)-H))) THEN 
AA = A(T2(2)-1) 
QA = Q(T2(02)-1) 
SA = S(T2(2)-1) 
ELSE 


22S EXTRAPOLATE TO LEFT FACE OF CONTACT SURFACE----- 


IF (12(2).EQ.2) THEN 
CALL BBORY(RR(1),QQ(1),S01),RRA,QQA,SA,AA,QA) 


ELSE 

CALL EXTRAP(RRIT2( 2 J-1),RRIT20 2 )-2),QQ(I2(02)-1),QQ(T2(2)-2), 
c S$(I202 J-1),S01202 )-2 ),XA02),H»RRA»QQA,SA,AA »QA ) 

END IF 

ENO IF 


SB = S(I2(2)4+1) 

CALL CSJUMP(AA,QA,SA,SB,G2,QB,AB) 
QQ@B = QQCALC(QB,AB,SB) 

RRB = RRCALC(QB,AB,SB ) 
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manne INTERPOLATE TO NODE(I2(2)) AND ASSIGN CORRECTED VALUES----- - 


aaa 


IF (12(2).EQ.N) THEN 
CALL BBDRY(RRB,QQB,SB,RRIT2(2)),QQ(1202)),S(1212));, 


C A(T2(2)),Q(T2(2))) 
ELSE 
CALL INTERP(RRB,RRI T2(2)+1) ,QQB,QQ(T2(2)+1),SB8> 
C S(T202)41),X8(2),(XB(2 +H) »,RRIT2(2)), 
Cc QQ(T2(2)),S(T202)),A0T202)),Q0T202))) 
END IF 
C 
C --=- SHOCK ON LEFT, HEADED RIGHT, JUMPS NODE WITH NO CONTACT --- 
C 
ELSE IF((SHOCK.EQ.321).AND.(CNTACT.EQ.100)) THEN 
I2(1) = NODE + 1 
CALL DELTAX(I2,SIGMA,XB,XA,H) 
C 
C e---- CHECK FOR A CONTACT SURFACE INTERACTING WITH SHOCK AT NODE----- 
c 
IF(C(SIGMA(2,2).LT.(X201)4H)). AND. 
Cc (SIGMA(2,2).GT.(X2(1)))) THEN 
AB = A(TI2(1))} 
QB = Q(I2(1)) 
SB = S$(T2(1)) 
ELSE 
C 
C eeee- EXTRAPOLATE TC RIGHT FACE OF SHOCK----- 
C 
IF (12(1).EQ.N) THEN 
CALL BBDRY(RRIN),QQ(IN),S(N),RRB,QQB,SB,AB,GB ) 
ELSE 
CALL EXTRAP(RR(I2(1)),RR(I2(1)41),QQ(1T2(1)),QQ(T2(1)41), 
Cc $(I2(1)),S(¢IZ2(1)+1),XB(1),H,RRB,QQB,SB,AB,QB) 
END IF 
END IF 
C 
C --- CALCULATE VARIABLE CHANGE ACROSS SHOCK --- 
G 
CALL SKJUMP( AB ,QB,SB,AR,DQ,VS>»G,G1»W,AA ,QA5SA ) 
QQA = QQCALC(QA,AA,SA) 
RRA = RRCALC(QA,AA,SA) 
C 
C e---- INTERPOLATE TO NODE(I2(1)-1) AND ASSIGN CORRECTED VALUES----- 
C 
IF (12t1).€Q@.2) THEN 
CALL BBDRY( RRA ,QQA,SA,RRIT2Z(1)-1),QQ(T2(1)-1)>, 
C S$(¢IT201 J-1),A0T 201 )-1),Q(IT201)-1)) 
ELSE 
CALL INTERP(RRA,RR(I2Z2(1)-2) ,QQA,QQ(I2(13-2),SA;, 
G SC I201)-2),XA(1),(0XA0C1I4H),RRII2Z(1)-1), 
C QQIT2(1L-1),S(I201)-1),A0CT2( 1-1); 
Cc Q(T2(1l)-1)) 
END IF 
cS 
C --- SHOCK ON RIGHT, HEADED LEFT, JUMPS NODE WITH NO CONTACT SURFACE-- 
C 
ELSE IF((SHOCK.EQ@.231).AND.(CNTACT.EQ.100)) THEN 
I2(1) = NODE 
CALL DELTAX(I2,SIGMA ,XB,XA,H ) 
IF( (SIGMA(2,2).GT.(X2(1)-(2.000#H))).AND.(SIGMA(2,2).LT. 
C (X201)-H))) THEN 
AA = A(T2(1)-1) 
QA = Q(I2(1)-1) 
SA = S(T2(1)-1) 
ELSE 
C 
C ----- EXTRAPOLATE TO LEFT FACE OF CONTACT SURFACE----- 
Cc 
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IF (132(1).EQ.2) THEN 
CALL EBORY(RR(1),QQ(1),S(1),RRA,QQA,SA,AA,QA) 
ELSE 
CALL EXTRAP(RR(I2(1)-1),RR(T2(1)-2) ,QQ(I2(1)-1),QQ(I2(1)-2), 
$(I2(1)-1),S(I201 )-2),XA(1),H> RRA, QQA,SA,AA>,QA ) 
END IF 
END IF 


w-- CALCULATE VARIABLE CHANGE OVER SHOCK --- 


IF(I2(1).EQ.N-1) THEN 


END IF 


CALL SKJUMP(AA,QA,SA,AR»,DQ,VS,G,G1,WN,AB ,QB,SB } 


IF(I2(1).EQ.N-1) THEN 


END IF 


END IF 


QQB 
RRB 


QQCALC(QB,AB,SB) 
RRCALC( QB ,AB »SB ) 


IF (I201).EQ.N) THEN 
CALL BBORY( RRB ,QGB,SB,RR(12(1)),QQ(I2(1)),S(T2(1)), 
A(T2(1)),Q(T2(1))) 

ELSE 

CALL INTERP(RRB,RR(I2(1)41),QQB ,QQ(I2(1)+1),S8, 
SCT201)41),XB(1),(X801)4+H) »RRIT2(1)), 
QQ(T201)),501T201)),A0CT2°01)),Q0T201))) 

END IF 


@-- CHECK IF SECOND NODE NEEDS CORRECTING IF SO LOOP BACK AROUND --- 


IF(RNODE(1).EQ.0) GO TO 40 


GO RETURN 
END 


NOOE = RNODE(1) 
SHOCK = RNODE( 2) 
CNTACT = RNOOE(3) 
RNODE(1) = 0 

GO TO 35 


SUBROUTINE DELTAX(1I2,SIGMA ,XB>XA>H) 


HE FEE IE IE IE IE IEE FETE IE IE IE IE IE IE IE IE IE JE FETE IE IE IETE JE IE IE IE FE IE IE IEDE JE IE JE JE IE IE IE IE IE IETE JE IE IE IE IEE JE IE IE 


x 


x 


* LOCATE DISCONTINUITY WITHIN INTERVAL SUBROUTINE ¥ 


> 


x 


FEM FEE HE IE IE IE IE IE HE IE IE IE IE IE HEHE IE IEE IE IE NE HE FE IE IE IE IE HE IE IE FE IE HE IE IE IE IE FE IE IE IE IE IE IE IE IE IE IE IE IEE IE HE 


DIMENSION 1I2(4) ,SIGMA(4,2 ),XB(@) ,XA(4),X2(4) 


INTEGER T2 

DOUBLE PRECISION SIGMA »>XB,XA>,H»X2 

X2(3 )=0.000 

%2(4 )=0.000 
——== CALCULATE DISTANCE NODE TO RIGHT OF DISCONTINUITY IS FROM LEFT -- 
<<2=ce BOUNDARY OF TUBE THEN DETERMINE DISTANCE FROM THIS NODE I2(*)---- 


----- TO DISCONTINUITY AND DISTANCE FROM NODE TO LEFT OF DISCONTINUITY- 
----- TO THAT DISCONTINUITY----- 


X2(1) 
*B(1) 
XACT) 


ee CALCULATE 


OBLE(T2(1)-1) * H 
(X201)-SIGMA(1,2)) 
(H-AXB(1)) 


SAME VALUES FOR CONTACT SURFACE----- 
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X2(2) = DBLE(I2(2)-1) ® H 

AB(2) = (X22 )-SIGMA( 2,2) 

XA(C2) = (H-XBi2)) 

RETURN 

ENO 

SUBROUTINE SKJUMP( ADN,Q0N,SON,ARATIO ,DELTAQ,VSHOCK »G,G1>H, 
Cc AUP ,QUP ,SUP ) 

388 BBB BUBB ERB BUBB BEEBE RBBB EBNEHE 

» ® 

xe SHOCK JUMP SUBROUTINE 

& * 


EIEN FETE IE IE HE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IEE IE IE TEE HE IE ETE HE IE HE PETES HE IEE HE 


*DN = VARIABLE DOWNSTREAM OF SHOCK 
ARATIO = SPEED OF SOUND RATIO 


D 


ELTAQ = VELOCITY JUMP ACROSS SHOCK 


*UP - VARIABLE UPSTREAM OF SHOCK 
VSHOCK = SHOCK VELOCITY 


DOUBLE PRECISION AON,QON,SON,ARATIO,DELTAQ,VSHOCK ,G,G1>W, 
Cc AUP ,QUP ,SUP ,SA1 ,SA2 ,UDN ,UUP 


AUP = ADNX¥ARATIO 


SA1 = (G1/G )*DLOG( (2. DOOXG#( W¥¥*2 )-G+1.000 )/1G+1.000)) 
SA2 = G1*DLOG(((G-1.000 }*(Wx*2 1423711641. D000 xt Wexe2 DD) 
SUP = SON - SAl - SA2 


UON = QDN - VSHOCK 

UUP = UON + ADNxXDELTAQ 
QUP = UUP + VSHOCK 
RETURN 

END 


SUBROUTINE CSJUMP(ADN,QON,SON »SUP ,G2 ,QUP ,AUP ) 


JE IE IE IE IE IE HE IE IE IE IE FE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IEIEIE JEHE JE IE IE IE IE IE TEI 


% % 
% CONTACT SURFACE JUMP SUBROUTINE ial 
x % 


JOBE EUR BEEBE BUBB RRHOHE 
DOUBLE PRECISION QUP ,QON,AUP ,ADN,SUP ,SON,G2 


CALCULATE THE SPEED OF SOUND CHANGE ACROSS THE CONTACT SURFACE--- 
NOTE: THE VELOCITY IS CONTINUOUS ACROSS A C.S.----- 


QUP = QON 

AUP = ADN * DEXP((SON-SUP )/G2) 

RETURN 

END 

SUBROUTINE EXTRAP(RRL >RRR>QQL ,QQR,SL>SR»,DX,DOH»RRINT »QQINT ,SINT >» 
Cc AINT »QINT ) 

JOB O UR BOEEIUE UBB O RRR EB BEB BRHEE 

~ * 

sal EXTRAPOLATION SUBROUTINE * 

* * 


IEE IEE IEE IEE HE IE IE IE IE IE IE IE IE IE IE IE IEIE IE JERE FE IE IEIEIE FE IE IEIE IE IE IE JEIE IE IE IE IE IE IE IE IE IE IE IE IE IEE TE IE IE 
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DOUBLE PRECISION RRL>RRR,QQL,QGR,SL,SR»,DX,OH »RRINT ,QQINT ,SINT 
c AINT ,QINT 


w----REIMAN VARIABLES ARE EXTRAPOLATED, AND THE CORRESPONDING----- 
woe--VALUES FOR VELOCITY AND SPEED OF SOUND ARE CALCULATED----- 


RRINT = RRL - (OX * (RRR=-RRL)/OH) 
QQINT = QQL - (DX * (QQR-QQL)/DH) 

SINT = SL - (0X * (SR-SL)/0H) 

AINT = (QQINT=-RRINT )/(2.000*SINT ) 

QINT = (QQINT+RRINT )/( 2.000) 
RETURN 
END 
SUBROUTINE INTERPC(RRL>RRR,QQL,QQR»,SL»SR»DX,D0H,RRINT ,QQINT ,SINT 5 

Cc AINT ,QINT } 

JE OE GOE IEE SEE IOC SIE FEE CRE IE RE EE FOIE IEE EE BEBE EEIOIE e  R 
x * 
* INTERPOLATION SUBROUTINE * 
x x 


HE FEE IEEE FE IE FE IE FETE IE IESE HE ESE SE ESE EE SE IESE SE 5E IE SE 36 3E IE IE IE IEE IE IE FETE IE IE IE IE IE IE IE IE JEIE IEIE IE IE HE 


DOUBLE PRECISION RRL>RRR>,QQL,QQR,SL>,SR»,DX,DH,RRINT ,QQINT ,SINT, 
Cc AINT ,QINT 


RRINT = RRL - (OX * (RRL-RRR)/DH) 
QQINT = QQL = (DX * (QQL-QQR)/DH) 


SINT = SL - (DX * (SL-SR)/DH) 
AINT = (QQINT-RRINT )/(2.D000*SINT } 
QINT = (QQINT+RRINT J/2.000 
RETURN 
END 
SUBROUTINE DBURST(N,H,QQ,RR»S,G,G61,62,DELT,I2,X2,W,AR,DQ,VS> 
RLWPRES ,SIGMA,A,Q) 
EH 3EIEIE 3E IE IE IE IE HE IE EI FE IE FE IE IEE FE IEIE HE IE FE IE TEI FE IE FE IEIE HE IE IE IE IE IE HE HEE HE IE IE NE HE 
% x 
€ DIAPHRAM BURSTING SUBROUTINE * 
% * 


FE IE IESE IE IEE HE IE IE IE IE TENE FE IE IE IEE FE FE IE IE HE IE IE IE IE IE IE FE IE IE IE IE IE IE IE IE IE IE IE IE IEE FE HEE 


INTEGER N>I,Y,I2,L»SHKDIR,CSDIR,LWPRES 

DIMENSION X2(4),RRON),QQ(N),S(N),I204) »SIGMA(G,2),A0N),Q(N) 

DOUBLE PRECISION X2,X%,H,AB,SA,SB,AA,QA ,QB ,QQA ,QQB,RRA,RRB,SIGMA, 
C RR »,QQ,S,TS»WN,DQ,AR»PR»G»G1,G2 ,VS,DELT »>CSRMN>A> 
Cc Q,MREIMN,DREMN,EREIMN WW SAL SA2 ,SAP »SBP > XB XA 


+++4++4+ LOCATING THE NODE TO THE RIGHT +++++4++4+++4++++ 


DO 10 L=1,;4% 
SIGMA(L»,1)=SIGMA(L>2) 
Y=0 
X=H 
T=2 
11 IF (.NOT.(Y.EQ.0)) GOTO 10 
IF (SIGMA(L,1).LT.X) THEN 
X20 L3=X 
T2(L)=I 
Y=1 
END IF 
X=X+H 
I=I+1 
GOTO 11 


10 CONTINUE ; 
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44444444+++++ CORRECT FOR DIAPHRAM BURSTING +4++44+4+44+44+4+4 


===-= AT TIME ZERO DETERMINE CORRECT SHOCK DIRECTION----- ~-- 
===-- SHKOIR = 3 IS A SHOCK HEADED LEFT, AND SHKDIR = 2 IS SHOCK----- -- 
= - === TRAVELING RIGHT----- ah. 


IF(LWPRES.EQ.3) THEN 


SHKDIR = 3 
X2(1) = X2(1) - H 
T2(1) = T2t1) - 1 
A2(2) = X202) - H 
IT2(2) 2 I2t2aee 1 
GO TO 20 

ELSE 
SHKOIR = 2 

END IF 


20 RRA=RR(T2(1)-1) 

RRB=RR(I201)) 
QQA=QQ(T2(1)-1) 
QQB=QQ(T2(1)) 
SASS(T201)-1) 
SB=S(T201)) 

21 AB=(QQB-RRB )/02.000%SB ) 
AA=(QQA-RRA J7(12.DOOXSA ) 
IF(SHKDIR.EQ.3) THEN 

MREIMN = (RRB-RRA)/AA 
DREMN = DABS(MREIMN) 
ELSE 
MREIMN = (QQA-QGQB )/AB 
DREMN = MREIMN 
END IF 


100 WHW=(3.0396408D01-( (DREMN+#2. 7574D00 )/0. 286337D00 ) ) 
W=5.513294000-DSQRT (HW) 
DQ=2.D00*( WeN-1.D00)/( W(G+1.D00) ) 
AR=DSQRT(2.000*%(G-1.D00 )*(1.D00+((G-1.000 )*WXW/2.D00 ) J* 

C (G¥G2*W¥W-1.D00) 70 (G41. D00 )¥W) 

PR=(2.D00*G/(G+1.D00 ) )*W¥W-( (G-1.000)/(G+1.D00)) 
OR=((G-1.000 )*¥W¥W4+2.D00 J/((G+1.D00 J) xWeW ) 
EREIMN=DQ+( AR-1.D00 )*G2-( ARXG1/G )¥DLOG( PR¥( DR¥XG ) ) 
IF (DABS(EREIMN-DABS(MREIMN)).LT.O.1D-5) GO TO 110 
DREMN = (DABS(MREIMN) - EREIMN) + OREMN 
GOTO 100 


110 SAl = (G1/G )¥DLOG( (2. DOOXG*( Wxx2 1-G+1.D00)/(G+1.D00)) 
SA2 = GI*XDLOG(((G-1.D00 )*(We*2 142)/70 (G41. D0O0 J¥( WR? ) ) ) 
IF(CSHKDIR.EQ.2) THEN 


SAP = SB - SAl =- SA2 
SBP = SAP 
ELSE 
SBP = SA - SAl = SA2 
SAP = SBP 
END IF 
103 IF(SHKDIR.EQ.2) THEN : 
CSRMN =((DEXP((SBP-SA)/G2 ) )*(SA )-(SBP ) JAR 
ELSE 


CSRMN =( ( DEXP((SAP-SB )/G2 ) }*(SB )-( SAP ) )*AR 
END IF 
IF(SHKDIR.EQ.3) THEN 
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MREIMN =( (RRB-RRA )/AA )+CSRMN 
DREMN = DABSIMREIMN) 


ELSE 
MREINN =( ( G@QA-QQB )/AB )-CSRMN 
DREMN = MREIMN ar 
END IF ome 
101 WAN=(3.0396408D01-( (DREMN42. 7574000 )/0. 286337000) ) 


W=5.513294D00-DSQRT (WH) 
DQ=2.D00*(WxH-1.000 )/( W(G+1.000)) 
AR=DSQRT(2.000%(G-1.D00 )¥(1.D004((G-1.D00 )*¥H*W/2 . 000 ) J 
Cc (G¥G2*WXN-1.D00) 70 (G+1.000 )*W) 
PR=(2.000*G/(G+1.000 ) )¥NxW-{ (G-1.000 )/(G+1.000)) 
OR=( (G-1.D00 )¥WN¥N+2.000 )/( (G41. D00 J¥WXW ) 
ERE IMN=0Q+( AR-1.000 }*G2-( AR*G1/G )x0LOG( PR¥( DRG ) } 
IF (DABS(EREIMN-DABS(MREIMN) ).LT.O.10-5)} GO TO 102 
OREMN = (OABS(MREIMN) = EREIMN) + OREMN 
GOTO 101 
102 SA1l = (G1/G )*DLOG( (2. DOOXG*(Wex*2 }-G+1.000 )/(G+1.000)) 
SA2 = G1lXDLOG( ((G-1. 000 )3( Wexe2 42 70 (G41. 000 D*( Wee?) )) 
IF(SHKOIR.EQ.2) THEN 
SAP = SB - SAl - SA2 
ELSE 
SBP = SA = SAl = SA2 

ENO IF 
IF (OABS(SAP-SBP).LT.0.10-5) GO TO 105 
IF(SHKDIR.EQ.2) THEN 

SBP = (SAP-SBP) + SBP 
ELSE 

SAP 
END IF 
GO TO 103 


(SBP-SAP)} + SAP 


DSSS CALCULATE SHOCK VELOCITY ----- 


105 IF(SHKOIR.EQ.3) THEN 
OQ = (-1.000)*DQ 
VS = CCRRA+QQA )*¥0.5D00) = (WAA) 
ELSE 
VS=(QQB+RRB )*0 . 5SDO0+H*AB 
END IF 
IF(SHKDIR.EQ.2) THEN 


= EXTRAPOLATE TO RIGHT FACE OF SHOCK----- 


CALL EXTRAP(RR(I2(1)),RROT2(1)41),QQ(T2(1)),QQ(IT2(1)+1), 
Cc S(I2(1)),S0I201)+1),H,H,RRB,QQB,SB,AB,GB ) 


CALL SKJUMP(AB,QB,SB,AR»,0Q,VS,G,G1 5W,AA ,QA,SA ) 


72-2 CALCULATE CHANGE IN VARIABLES ACROSS CONTACT SURFACE----- 


QB = QA 
AB = AA 
SB = SA 
SA = S(T202)-2) 


‘CALL CSJUMP( AB ,QB,SB,SA,G2,QA,AA ) 
QQA = QA + AAXSA 

RRA = QA - AAX®SA 

XA=0.00 


CALL INTERP(RRA,RRIT2( 2 1-2) ,QQA,QQ(T2(2)-2),SA5 
S(¢T2(2 )-2) XA, (XA4H) RRO IT2( 2 )-1)5 
QQ(T2(2)-1),S0T202 }-1),ACT202)-1), 
Q(T202)-1)) 


oan 


ELSE 
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weem-EXTRAPOLATE TO LEFT FACE OF SHOCK----- 


Cc 


=--=-CALC 


CALL EXTRAP(RR(I2(1)-2),RRII2(1 )-2),QQ(I2(1 )-1),QQ(I2(1)-2), 
$(I201)-1),S(12(1)-2),H»H, RRA, QQA,SA,AA,QA) 


ULATE CHANGE IN VARIABLES ACROSS SHOCK----- 


CALL SKJUMP(AA,QA,SA,AR ,0Q,VS,G,G1,W,AB,QB>SB ) 


w----CALCULATE CHANGE IN VARIABLES ACROSS CONTACT SURFACE----- 


QA = QB 

AA = AB 

SA = SB 

SB = S(T2(2)+1) 

CALL CSJUMP(AA,QA,SA,SB,G2,QB,AB ) 
QQ@B = QB + AB*SB 

RRB = QB - AB*SB 

XB=0.00 


Re ame INTERPOLATE TO NODE(I2(2)) AND ASSIGN CORRECTED VALUES----- 
CALL INTERP(RRB,RR( 1202 }+1),QQB,QQ(I2(2)4+1),SB8, 
C S( 1202 )41),XB,(XB+H) »,RRIIT202)), 
C QQIT202)),S01202)),A0T202)),Q0T202))) 
END IF 
RETURN 
END 
SUBROUTINE BONDRY(Q1,Q2,QB,A1,A2,QQ1,QQ25,RR1»RRZ,S1,S2>H>EE, 
C DELT,BNORY ,BOPRS ,BOPR,BDDR,BOTR>J> 
C NEWQQ1 ,NEWRR1,NEWS1»,G,G1,G2,HALT ,BORY SK ) 
FOIE IE IE IE IE IE IE IE IE IE SEI IE IE 3 IE 36 I IE IE EE EE IE IE 336 36 36 8 IE IE 3 I 38 36 3 96 38 36 3b IEE IE I 38 IE 6 36 3 36 38 EE 
*% % 
% BOUNDARY COUNDITION SUBROUTINE % 
*% % 
JERE EEE EERE BEBEUHBHEHE 
2s == === S=-— VARIABLE DEFINITIONS ------------ 
*®BD - VARIABLE AT PHANTOM NODE 
%*1 - VARIABLE AT BOUNDARY NODE 
*%2 =- VARIABLE AT FIRST NODE INSIDE BOUNDARY 
D2 = DENSITY RATIO AT FIRST NODE INSIDE BOUNDARY 
P2 = PRESSURE RATIO AT FIRST NODE INSIDE BOUNDARY 
QB - INITIAL VELOCITY AT AN OPEN BOUNDARY 
T2 = TEMPERATURE RATIO AT FIRST NODE INSIDE BOUNDARY 
DIMENSION AINT(3),2(3),QPRIM(3),APRIM(3 ),INTEG(3),AAVG(3), 
QQBD(6),RRBD(6),SB0(6),QBD0(6),ABD( 6) 


INTEGER BNDORY ,BOPRS>J,HALT »SK,BORY »K,L»M>T 


DOUBLE PRECISION RRBD,QQBD,SB0,QBD ,ABD »BOPR,BODR »BOTR »DELQQH , 


Oaaaan 


DELQQL ,DELRRH,DELSH,DELSL,DELQH,DELAH,DELQL, 
DELAL ,»H,EE,DELT »,QQINT ,RRINT »>SINT »QPRIM,APRIM, 
AINT 52 »SAVG, AAVG ,Q1 ,Q2Z,A2,A1L,RR1»,QQ15S1,QQZ2,RR2, 
S2,D0LTAQQ,OLTARR ,»,OLTAS,INTEG,QQSTEP »RRSTEP,D2;, 
SSTEP »,NENQQ1 »NEWRR1,NEWS1,G,G1,G2,DELRRL,T2,P2, 
AR ,O0Q,VS,W,QB,NENTR »NENDR,NEWPR ,EEL,EEZ,QRR»QQQ 


COMMON AR,0Q,VS,W 


----- DET 


ERMINE CURRENT PRESSURE AT NODE JUST PRIOR TO BOUNDARY ----- 


((QQ2-RR2 )**2 371%. DOOKI S22 ) 3} 
((1.000/T2 J¥DEXP(Ge(1.D000-G )*(S2-G2 ) ) )¥*(-G1) 
T2 * D2 
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C ----- DETERMINE IF LEFT OR RIGHT BOUNDARY ----- 
Cc 
IF(BORY.EQ.%) THEN 


ELSE 


zrexXK ee ls 


£ Ul oo Wn 


END IF 
IF((J.EQ.1).AND.(BNDRY.EQ@.0).AND.(BOPRS.EQ.0)) THEN 
QBD(K) = 0.000 

ENO IF 

QQBO(L) 
QQBD(M) 
RRBO(L ) 
RRBO(M) 
SBD(L) 

SBO(M) 

QBD(L) 

QBD(M) 

ABDO(L) 

ABD(M) 


ao ot 
a 
a 
ps 


ob 0 te kt 
&) 
p<! 


a —— DETERMINE IF BOUNDARY IS OPEN OR CLOSED, AND SET PROPER ----- 
C 2---- VALUES AT PHANTOM NODE “LBD”. BNORY = 1(CLOSED),O( OPEN) ----- 


5 IF (BNORY.EQ.1) THEN 
RRBO(K) = -QQBO(M) 
DABS( RRBO(M)) 


W 

w 

oO 

~ 

om 
woo OF 

W 

@ 

oO 

= 


ELSE IF (BOPRS.EQ.1) THEN 
IF (J.EQ.1) THEN 
RRBO(K) = RRBDIL) 
QQBD(K ) QQBO(L) 
SBO(K) SBDO(L) 
ABO(K ) ABD(L) 
QBD(K ) QB 


END IF 
GO TO 10 
ELSE IF ((BOPR-P2).LE.0.1000) THEN 
Cc IF(SK.EQ.3) GO TO 35 
SBD(K) = SBD{L) 
BOOR = ((1.0DD0/B0PR )¥( DEXP( (SBD(K )-G2 )*(-(G/G1))))) 
c **(-(1.000/G)) 
BOTR = ((BODR )**{1.000/G1) )*DEXP(G¥*(1.000-G )*( SBD(K )-G2)) 
RRBDIK) = QBO(K )-( DSQRT(BOTR ) )*SBO(K ) 
QQBD(K) = QBD(K)+(DSQRT(BOTR) J*SBD(K) 
ABD(K) = (QQBD(K) - RRBD(K)) / (2.D00*SBDIK)) 
RRBD(L) = RRBO(K) 
QQBD(L) = QQBDIK) 
SBD(L) = SBO(K) 
QBO(L) = QBO(K) 
ABD(L) = ABD(K) 
GO TO 10 


ELSE 
PRINT % »*'THE OPEN BOUNDARY HAS A PRESSURE HELD CONSTANT ' 
PRINT * »'THAT IS HIGHER THAN THE PRESSURE INSIDE THE TUBE' 
PRINT * ,"ADDITIONAL CODE IS NEEDED TO PRECEDE' 

, HALT = 1 
GO TO 20 

ENO IF 


C e---=- CALCULATE THE JUMP TO THE NEXT TIME STEP USING CHARACTERISTICS-- 
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Cc 


C --- CHECK FOR SHOCK LEAVING TUBE OR WITHIN H OF BOUNDARY AND PICK -- 
C --- CHOSE APPROPRIATE ALGORITHM TO ADVANCE BOUNDARY TO NEXT TIME --- 
c 


10 IF(SK.EQ.2) THEN 
IF(BORY.EQ.2) THEN 
CALL SKJUMP(ABD(L41),QB0(L4+1),SB80(L41),AR,0Q,VS,G6,G15N, 
Cc ABO(L),QBD0(L),SBD(L)) 
QQBD(L) = QBO(L) + ABO(L)*SBDO(L) 
RRBO(L) = QBO(L) - ABD(L)*SBDO(L) 
GO TO 35 
ELSE IF(BORY.EQ.3) THEN 
CALL SKJUMP( ABD(L-1),QBD(L-1),SBD0(L-1),AR,0Q,VS,»,G,G1>H, 
C ABD(L),QBD(L),SBO(L)) 
QQBD(L) = QBD(L) + ABD(L)*SBDO(L) 
RRBO(L) = QBD(L) - ABDC L)¥SBDO(L) 


GO TO 35 
END IF 
ENO IF 
DELQQL = QQBD(L) - QQBD(L-1) 
DELQQH = QQBD(L+1) - QQBDI(L) 
DELRRH = RRBD(L+1) = RRBO(L) 
DELRRL = RRBO(L) - RRBO(L=-1) 
DELSH = SBD(L+1) - SBD(L) 
DELSL = SBD(L) - SBD(L-1) 
DELQH = QBD(L4+1) - QBD(L) 
DELQL = Q@BD(L) - QBO(L-1) 
DELAH = ABD(L4+1) - ABD(L) 
DELAL = ABD(L) - ABDIL-1) 


IF( (BNORY.EQ.0}.AND.(SK.EQ.1)) THEN 
IF (BORY.EQ.2) THEN 
CALL COND3(QB0(L?),QBO(L4+1),ABD(L),ABD(L4+1),RRBO(L),QQBD(L), 


Cc SB8D(L),DELQQH,DELRRH,DELSH,DELQH,DELAH ,DELT »H>EE > 
Cc QQINT ,RRINT »>SINT ,QPRIM,APRIM,AINT ) 

ELSE 

CALL COND2( QBD0(L),QBO(L-1),ABD(L),ABO(L-1),RRBO(L),QQBD(L), 
Cc SB80(L),OELQQL,DELRRL,DELSL,DELQL,DELAL ,DELT,H,EE, 
Cc QQINT »>RRINT »>SINT »>QPRIM »APRIM >AINT ) 

END IF 


ELSE IF((8NORY.EQ.1).AND.(SK.EQ.1)) THEN 
QQSTEP = 0.000 
RRSTEP = 0.000 
SSTEP = 
GO TO 30 


oS ----— USE CONDITION 1 ALGORITHM TO CALCULATE RRINT ,QQINT>,SINT ----- 


ELSE 

CALL COND1( QBD0(L),QBO(L+1),AB0(L),ABD(L4+1),RRBO(L),QQBD(L); 
C SBD(L}),DELQQL,DELRRH,DELSH,DELSL,> 
Cc DELQH »DELAH,DELQL ,DELAL,H,EE ,DELT ,QQINT »RRINT »SINT > 
Cc QPRIM>,APRIM >AINT ,QBD(L-1 ),ABD(L=-1)) 


eevcncece CALCULATE OLTA QQ, DLTA RR & DLTA S 


aan 


DLTAQQ=QQINT-QQBDI(L) 
DLTARR=RRINT-RRBD(L) 
DLTAS=SINT-SBO(L) 


------ CALCULATE 2(K)'S ----------- 


aaan 


AAVG(1 )=(AINT(1)4AB0(L ))72.0000 

AAVG(3 )=( AINT(3 )+ABD(L))72.0000 

AAVG( 2 )=0.0000 

SAVG = (SINT+SBO0(L))72.000 
Z(1)=-(1.0000/G2 )¥AAVG(1 )*( SAVG-G2 )*( QPRIM( 1 )-G2*APRIM(1 )) 
Z2(3 )=(1.0000/G2 JXAAVG(3 )*( SAVG-G2 )¥(QPRIM( 3 )+G2*APRIM( 3 )) 
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2(2)=0.0D00 
----- INTEGRATE THE Z(K)'S ------------ 


INTEG( 2)=0.0D000 
INTEG(1 )=2(1 )*XOELT 
INTEG(3 )=Z2(03 J*DELT 


w--2---- SOLVE THE EQUATION ---------- - 


QQSTEP=DLTAQQ+INTEG( 1) 
RRSTEP=DLTARR+tINTEG(3 ) 
SSTEP=DLTAS+INTEG( 2) 


30 NEWNQQ1=QQBD(L )+QQSTEP 
NEWRR1=RRBD(L)+RRSTEP 
NEWS1=S8D(L)+SSTEP 


@--- UPDATE PHANTOM NODES FOR OPEN BOUNDARY CONDITION ----- 


35 IF ((8NDRY.EQ.0).AND.(8DPRS.EQ.1) THEN 
RRBD(K) = RRBD(L) 
QQBDIK) = QQBDdD(L) 


SBD(K) = SBDIL) 
ABD(K) = ABD(L) 
QBD(K) = QBD(L) 


END IF 
IF ((8NDRY.EQ.0).AND.(BDPRS.EQ.0)) THEN 

IF (SK.EQ.2) THEN 

QBD(K) = (QQBDILI+RRBD(L)) 7 2.000 

ELSE 

NEWTR= ((NEWQQ1 - NEWRR1 )**2)/(04.D00%( (NEWS1 )¥*2 ) ) 

NEWDR= ((1.D00/NEWTR )¥DEXP( G¥(1.D00-G )*(NEWS1-G2 ))) 

Cc #¥(-G1) 

NEWPR= NEWTR#¥NEWDR 

EE1 = DA8S( BDPR-NEWPR ) 

EE2 = DA8S(SB8D(L)-NEWS1 ) 

IF ((EE1.GT.0.1D-5).OR.(EE2.GT.0.10-5)) THEN 
QRR= NEWRR1I + (DSQRT(BDTR ) I*SBD(L) 
QQQ= NEWNQQ] - (DSQRT(BDTR) )*SBDIL) 
QBD(K) = (QRR+QQQ )/2.D0DD 
GO TO 5 

END IF 

END IF 

END IF 
20 QQl = QQBd(L) 
RR1 = RRBD(L) 


S1 = SBDIL) 

Ql = QBDIL) 

Al = ABD(L) 

RETURN 

END 

SUBROUTINE SRFLCT(QQA,RRA,SA,SIGMA,VS,DELT »LWPRES >RRB »QQB SB QB, 
Cc A8,G,G61,G2) 

38888008 RR BEEBE EEUU BUREN BUBBREHE 

* * 

* SHOCK REFLECTION AT SOLID BOUNDARY SUBROUTINE * 

x x 


HE HE IE SEIE HE FE IE IE IE IE EIE IE FE IE IE IE FE IE IE IE FE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE FE IE IE IE TE IE TE IE IE IE IE ETE IEE TE HEE 


------------ VARIABLE DEFINITIONS ------------ 


DELTEX - THE EXCESS TIME IN A TIME STEP WHEN THE SHOCK IS 
EXACTLY AT THE SOLID WALL 
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DELTWL =- THE TIME FOR THE SHOCK TO REACH THE WALL 


DIMENSION SIGMA(4 52) 

INTEGER LWPRES 

DOUBLE PRECISION QA,Q8,AA,AB,SA,SB,QQA ,QQB,RRA,RRB,SIGMA, 
C DQ,W,AR,PR,DR,EREIMN,SA1,SA2,VS,DELT ,DELTWL, 
Cc DELTEX,G»G1,G2 


2o--- CALCULATE THE TIME FOR SHOCK TO REACH WALL, AND EXCESS TIME ---- 
vores IN THIS TIME STEP ----- 


IF (LWPRES.EQ.2) THEN 
DELTWL = (1.000-SIGMA(1,1))/VS 
DELTEX = DELT - DELTWL 


cco-- CALCULATE VELOCITY GRADIENT OVER SHOCK AS IT REFLECTS ----- 


QA = (QQA+RRA)/2.D00 

QB = 0.000 

AA = (QQA-RRA)/(2.D00*SA ) 
DQ = (QB-QA)/AA 


W = DSQRT(( (DQ**2 )*0.36D00)+1.000) - (D0Q*0.6D00) 
we--- FOR SHOCK AT LEFT BOUNDRY DO THE SAME ----- 


EESE 
QQB=QQA 
RRB=RRA 
SB=SA 
DELTWL = (DABS(SIGMA(1,1))-0.000)/DABS(VS ) 
DELTEX = DELT - DELTWL 


———— CALCULATE VELOCITY GRADIENT OVER SHOCK AS IT REFLECTS ----- 


QB = (QQB+RRB)/2.000 

QA = 0.000 

AB = (QQB-RRB)/(2.D00%*SB ) 
DQ = (QA-QB )/AB 


----- CALCULATE EXACT REIMAN VARIABLE JUMP FROM DQ ----- 


W = DSQRT(( (0Q**2 )*0.36000)+1.000) + (0Q*0.6000) 
ENO IF 


i CALCULATE AR,PR,DR OVER SHOCK ----- 


AR=O0SQRT(2.000*(G-1.D00 J)*(1.000+((G-1.000 JxWNxW/2.D00 ) )® 
C (G*G2*W*W-1.000))/((G+1.000 JXW) 

PR=(2.000*G/(G+1.000 ) )x¥WxW-(1G-1.000 )/(G+1.000)) 

OR=((G-1.000 )¥W*XN+2.000 J/( (G+1.D00 JxWXW ) 

SAl (G1/G }*OLOG( ( 2. D00*G*( W*x2 )-G+1.000 )/(G+1.000)) 

SA2 G1*¥DLOG( ((G-1. DOO )*( Wx*2 )42 1/0 (G41. 000 XC Wx*2 ))) 


IF (LWPRES.EQ.2) THEN 
SB = SA - SAl = SA2 
AB = AA*AR 


-—-—= CORRECT REIMAN VARIABLES AT BOUNDARY ----- 


RRB 
QQB 


QB - AB*SB 
QB + AB*SB 


ro ga CALCULATE NEW SHOCK SPEED ----- 
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VS = ((RRA+QQA )*0.5000) - (HAA) 
----- CALCULATE POSITION OF REFLECTED SHOCK AT END OF THIS TIME STEP-- 


SIGMA(1,2) = VS*DELTEX + 1.000 


LWPRES = 3 
SSS SHOCK REFLECTING AT LEFT BOUNDARY ----- 
ELSE 
SA © SB - SAl =- SA2 
AA = ABXAR 
sccce CORRECT REIMAN VARIABLES AT BOUNDARY ----- 


RRA = QA - AA*SA 
QQA = QA + AAXSA 


ewere> CALCULATE NEW SHOCK SPEED ----- 
VS = ((CRRB+QQB )¥0.5D00) + (W¥AB) 
ew---- CALCULATE POSITION OF REFLECTED SHOCK AT END OF THIS TIME STEP-- 


SIGMA(1,2) = VS*DELTEX 
LWPRES = 2 
RRB=RRA 
QQB=QQA 
SB=SA 
QB=QA 
AB=AA 
END IF 
RETURN 
END 


SUBROUTINE BBDRY(RRI ,QQI1»,S1 »RR2»,QQZ S52 ,A2,Q2Z ) 


JEIEIETE FE IE IE IE JERE IE IE IE IE IE IE IE IE IE IE FETE IE IE JE IE 3G IE IE IE IE IEIE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE 


* x 
% NODES NEAR BOUNDARY SUBROUTINE x 
* * 


HE IEIEIE IE FE IE FETE FE IE IE IE IE JEIE IE FE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE HE IE IE TE IE 


DOUBLE PRECISION RR1,QQ1,51,RR2,QQ2,52 ,A2 QZ 


RR2 = RR1 
QQz = Qaal 
S2 = S31 


A2 = (QQ2 = RRZ2)7(2.D00%S2 ) 
= (QQ2 + RR2)/72.000 


SUBROUTINE BORDER( JSTOP ) 


FOIE FEE FE IE IE IE HE IE DE FE HE IE DE IE IE IE IE IE IE IE IE AE 3G IE IE IE IE IE IE IE IE IE IE IE IE IE IE FETE HE IE IE IE IE IE IE 


* * 
% PLOTTING AREA SETUP ROUTINE # 
* % 


HIE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE IE FE IE IE IE IE IE IE IE IE IE IE HE 


INTEGER I 

REAL XORG(4)/1.7554.6551.75,4.65/ 

REAL YORG(4)/2.75,2.75,5.80,5.80/ 

REAL YMAX(4)/5.5,5.5,1.5,6.20/ 

REAL YMIN(4)/0.55,0.55-0.554%.90/ 

DO 70 I=1;4 as 
CALL PHYSOR( XORG(I),YORG(I)) 
CALL NOBROR . 
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CALL AREA2D(2.4 2.4) 
CALL FRAME 
CALL GRAF(O.,'SCALE’,1.0,YMIN(I), ‘SCALE* ,YMAX(T)) 
CALL ENDGRIO) 
70 CONTINUE 
CALL PHYSCR(1.25,2.25) 
CALL NOBROR 
CALL AREA2D(6.00,6.50) 
CALL HEADIN( ‘SHOCK TUBE RESULTS$',100,1.7,4) 
CALL HEADIN( 'FIRST ORDER N =101$°,100,1.254¢) 
CALL HEADIN( ‘DENSITY RATIO = 5.0 TEMP RATIO = 
CALL HEADIN( ‘PRESSURE RATIO = 5.0$',100,1.0,4¢) 
CALL GRAF(O.,‘'SCALE',1.,0., 'SCALE' ,JSTOP } 
CALL ENOGRI(O) 
RETURN 
ENO 


1$',100,1.054) 


an 


SUBROUTINE PLOT(J>,JSTOP»N»,QQ,RR»S»H»XARRAY y 
#P ARRAY »DARRAY »QARRAY »SARRAY »G »G1»G2 ) 


3 36 36 36 3 E30 FE JE EE 38 EEE EE EE 3638 FEE 98 3058 390 38 06 9 38 96 9 8 EE EE 
* * 
* GRAPHICAL PLOTTING ROUTINE * 
* * 
20RD BORER BBE BBR HBBE 


NOAAaHAaAaANH 


INTEGER I,»N»,J,JSTOP,KNT(4 3/1,4,6,97 
DIMENSION QQ(N),RRIN)»SON),XARRAY(N) ,QARRAY(N) ,PARRAY(N), 
Cc DARRAY(N),SARRAY(N ) 
DOUBLE PRECISION QQ>RR>S+»H»G>,G1>,62 
REAL XARRAY ,QARRAY ,PARRAY ,DARRAY 
REAL SARRAY,TEMP,HR,GIR,»,GR»,G2R 
REAL XORG(4)/1.75,%4.6551.7554.65/ 
REAL YORG(4)/2.75,52.75,5.80,5.80/ 
REAL YMAX(4)/5.5,5.5,1.0,6.20/ 
REAL YMIN(G )/0.5,0.559-1.054. 907 
CHARACTER¥®@ IYNAM 
DIMENSION IYNAM(13 ) 
DATA IYNAM/ 'PRES‘, ‘SURE','$ ",'DENS', 
w°LUYS 5 “VEUO"°S GLiy ss *,'MODOI','FIED',* ENT’, ‘ROPY','$ FA 
Cc 
C --- CONVERT DOUBLE PRECISION TO SINGLE PRECISION <--- 
8 
GR=SNGL(G) 
GIR=SNGL(G1 ) 
G2R=SNGL(G2 ) 
HR=SNGL(H) 
DO 80 I=1>N 
QARRAY (TI J=(SNGL(QQ(T J4RRIT31/2.03 
TEMP=SNGL(QQ(T )-RR(TI ) J¥SNGL( QQ(T J-RR(I))/(4. OXSNGL(SCT )¥S(T)) 3 
OARRAY( TI )=((1.0/TEMP )¥EXP( GR*(1.0-GR )*(SNGL(S( I) )-G2R ) ) )¥¥(-G1R) 
SARRAY(T J=SNGL(S(T)) 
PARRAY( TI J=TEMP*#DARRAY(T ) 
80 CONTINUE 
0O 83 I=15,4¢ 
CALL PHYSOR( XORG(I),YORG(I)) 
CALL AREA2D(2.4,2.4¢) 
CALL XNAME( 'X',1) 
CALL YNAME( IYNAM(KNT(I)},100) 
CALL GRAF(Q.,‘SCALE',1.0,YMIN(I)»,°SCALE * »YMAX(T)) 
IF (I.EQ.1} CALL CURVE( XARRAY ,PARRAY ,N,0) 
IF (I.EQ.2) CALL CURVE( XARRAY ,DARRAY »N,0 ) 
IF (I.EQ.3) CALL CURVE( XARRAY ,QARRAY ,N,O ) 
IF (I.EQ@.4) CALL CURVE( XARRAY ,SARRAY,N,0O} 
CALL ENDGR(O} 
83 CONTINUE he 
RETURN 
ENO 
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aAaAnNgaAnaAN 


SUBROUTINE EXACT(N,.<INIT,7T,VHEAD »VTAIL,VCDE ,VSE ,OLCD,0LSH,QQ,RR;S, 
#H »>XARRAY ,OARRAY »G,G1,G2,ORI ) 


FE HEIEIE FETE HE IEIE IE HE JEIE JEIE IE JE IE IEE ETE DE IE IE FE IE IE IE I IEAE TE IE IE EIEIE HE IE IE IE IEIE TEE ETE 


% % 
a EXACT SOLUTION COMPARISON ROUTINE * 
x % 


FEI IEE EIEN IEE IE IE HE DE IE IE JE IE IE IE IE IEE IE IE HE HEIE IE IE HE IEIEIE HE ETE JE IEIE IE IE IE IE IE JEIE TIE TE 


INTEGER N 

DIMENSION Q@Q(N),RR(N},SON) ,XEXACT(6),YEXACT(6),XARRAY(N)>» 
c DARRAY(N ) 

DOUBLE PRECISION QQ,RR>S>H>G>G1,G2,T,DRI 

REAL XEXACT > YEXACT »XAINIT »VHEAD,VTAIL ,VCDE 

REAL VSE ,DLCO,OLSH,TEMP ,XARRAY ,DARRAY 


XEXACT(1)=0.0 

XEXACT ( 2 J=XINIT4+THVHEAD 
XEXACT(3 J=XINIT+T*¥VTAIL 
XEXACT (4 J=XINIT+T*¥VCDE 
REXACT(S JSEXINIT4+THVSE 
XEXACT(6 }=1.0 
YEXACT(1)=SNGL(ORT } 
YEXACT( 2 J=YEXACT(1 } 
YEXACT(3 )=OLCO 

YEXACT(4 J=YEXACT(3) 
YEXACT(5 )=DLSH 

YEXACT(6 }=1.0 


DO 90 T=1,5N 

QQ(T J=SNGLIQQ(T)) 

RR(UI J=SNGLORR(TI 7) 

S(T JSSNGLES(TI)) 

TEMP=(QQ(T )-RRUI ) *(QQ0TI-RRII) 704%. 0%S(TI *S(T )) 

DARRAY(T J=((1.0/TEMP )*EXP( G*(1.0-G )¥(SCI )-G2 )) )¥*(-G1) 
90 CONTINUE 

CALL PHYSOR(2.65,2.25) 

CALL NOBROR 

CALL AREA2D(6.75,4.0) 

CALL XNAME('X',1) 

CALL YNAME( 'DENSITY',7) 

CALL HEADIN( 'DENSITY OISTRIBUTIONS' ,100,1.3,4) 

CALL HEADIN( 'FIRST ORDER N =1015* ,100,1.0,4) 

CALL HEADIN( 'DENSITY RATIO = § TEMP RATIO = 1$',100,1.054) 

CALL HEADIN( ‘PRESSURE RATIO = 5$',100,1.0,4) 

CALL LINES( "EULER-1 SOLUTIONS ' ,IPKRAY »2) 

CALL LINES( "EXACT SOLUTIONS' ,IPKRAY,1) 

CALL FRAME 

CALL GRAF(0.0,'SCALE',1.0,0.0,'SCALE',6.0) 

CALL MARKER(O) 

CALL CURVE( XEXACT ,YEXACT »6,-1) 

CALL LEGLIN 

CALL CURVE( XARRAY ,DARRAY »N,0 ) 

CALL LEGEND(IPKRAY »,2,4.25,2.75) 

CALL ENDGR(O} 

RETURN 

END 
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